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ABSTRACT 


Current  limitations  of  atmospheric  predictive  skill  are  investigated  through  comparison 
of  correlation  and  error  statistics  of  operational  and  research  global  models  for  two  boreal  win¬ 
ter  seasons.  In  1993,  bias-corrected  models  produced  anomaly  correlations  of  0.6  after  7  days, 
with  relatively  little  forecast  skill  beyond  that  point.  In  2003,  the  forecast  skill  of  a  more  devel¬ 
oped,  higher-resolution  operational  model  has  been  extended  36-48  h,  while  the  skill  of  the 
unchanged,  low-resolution  research  model  has  been  extended  6-16  h.  This  implies  more  pre¬ 
dictable  patterns  in  2003  or  model  and  initial  state  improvements  made  since  1993.  The  rela¬ 
tive  importance  of  improved  model  resolution/physics  and  initial  state  to  the  lengthening  of 
forecast  skill  is  diagnosed  through  the  evaluation  of  root  mean  square  (rms)  evolution  of  ana¬ 
lyzed  and  forecast  differences  of  500-mb  height  and  meridional  wind.  Results  indicate  that 
forecast  sensitivity  to  initial  data  is  less  important  than  is  the  sensitivity  to  the  model  used. 
However,  the  sensitivity  to  model  used  (rms  of  model  forecast  differences)  is  comparable  or 
smaller  than  the  rms  forecast  error  of  the  model,  indicating  model  forecasts  are  more  similar  to 
each  other  than  to  reality.  Correlations  of  model  errors  to  each  other  quantify  this  similarity, 
with  correlations  exceeding  the  asymptotic  value  of  0.5  through  the  14-day  forecasts.  Investi¬ 
gations  of  initial  state  error  evolution  by  wavenumber  show  long  waves  generally  account  for 
more  of  the  total  uncertainty  growth  in  14-day  research  model  integrations  than  do  short 
waves,  with  the  relative  importance  of  each  determined  by  initial  state  uncertainty  and  model 
resolution.  Results  indicate  current  predictive  skill  may  be  impacted  by  model  sophistication, 


but  error  pattern  similarities  suggest  a  common  deficiency  of  models,  perhaps  in  the  initial 
state  uncertainty. 
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CHAPTER  1 


INTRODUCTION 

A  number  of  atmospheric  forecast  models  have  been  developed  and  now  display  con¬ 
siderable  skill  in  weather  prediction.  The  underlying  philosophy  of  these  developments  is  that 
improved  models  and  more  accurate  initial  conditions  should  provide  better^ forecasts.  The 
importance  of  specific  model  improvements  relative  to  specific  observational  enhancements 
may,  nonetheless,  still  be  inadequately  understood.  White  et  al.  (1999)  addressed  some  of  these 
questions,  and  suggested  the  error  spread  among  models  of  very  different  configuration  and 
resolution  is  generally  less  than  the  magnitude  of  the  error  in  any  single  relatively  advanced 
model.  Increased  resolution  in  regional  models  does  lead  to  some  improvement  in  skill,  partic¬ 
ularly  for  heavier  categories  of  precipitation  that  are  not  simulated  at  coarser  resolution  (White 
et  al.  1999),  and  recent  studies  show  improvement  in  wind  forecasts  due  to  enhanced  model 
resolution  (Hart  2004).  This  benefit,  however,  appears  to  be  relatively  small  for  other  forecast 
variables,  and  is  evident  only  for  short-term  forecasts  in  which  point  validations  are  made 
against  observations,  even  in  the  most  highly  resolved  and  dynamically  sophisticated 
approaches. 

Various  hypotheses  have  been  proposed  to  explain  this  result.  Some  of  these  note  cur¬ 
rent  observing  systems  contain  inadequate  resolution  of  local  regions  of  pronounced  dynamic 
instability.  Other  explanations  of  the  marginally  superior  performance  of  high-resolution  mod¬ 
els  point  to  the  evident  difficulties  in  the  validation  of  forecast  features  against  coarsely  spaced 
observations.  All  but  one  of  the  models  studied  by  White  et  al.  (1999)  were  limited  area  mod- 
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els  and  these  are  known  to  be  strongly  influenced  by  lateral  boundary  conditions  supplied  at 
their  perimeter  (Warner  et  al.  1997;  Paegle  et  al.  1997).  The  Medium  Range  Forecast  (MRF) 
model  was  the  only  global  model  evaluated  by  White  et  al.  (1999)  and  that  model  also  pro¬ 
vided,  or  strongly  influenced,  the  lateral  boundary  conditions  of  the  tested  limited  area  models. 

The  initial  stages  of  the  present  research  compared  two  global  models  [the  Utah  Global 
Model  (UGM)  and  the  MRF],  and  followed  a  hypothesis  suggested  by  Miguez-Macho  and 
Paegle  (2000,  henceforth  MMP).  This  perspective,  which  is  based  upon  downscale  uncertainty 
growth,  is  supported  by  MMP  and  by  Miguez-Macho  and  Paegle  (2001),  and  is  rooted  in  the 
early  barotropic  model  predictability  studies  of  Thompson  (1957)  and  Lorenz  (1969).  Results 
obtained  by  MMP  contrast  with  other  recent  literature  emphasizing  the  importance  of  rela¬ 
tively  smaller-scale  instabilities  (e.g..  Palmer  et  al.  1998;  Hartmann  et  al.  1995)  of  the  initial 
state,  and  other  local  error  sources  (e.g.,  Rabier  et  al.  1996).  These  perspectives  and  contrasting 
theories  are  summarized  by  MMP,  who  presented  preliminary  work  suggesting  that  the  domi¬ 
nant  source  of  short-term  forecast  errors  may  be  the  uncertainty  of  relatively  large  scales  of  the 
initial  state. 

The  present  research  extends  MMP’s  studies  to  sort  out  the  relative  roles  of  the  model¬ 
ing  technique  and  initial  state  uncertainty  and  continues  to  probe  the  limitation  to  deterministic 
weather  prediction  due  to  inadequate  observation  of  relatively  large  scales  of  the  atmosphere. 
Other  global  model  studies  attempting  to  prioritize  the  relative  contributions  of  initial  errors, 
boundary  errors,  and  model  errors  to  total  forecast  error  include  work  by  Reynolds  et  al.  (1994) 
and  Hacker  et  al.  (2003).  Simmons  and  Hollingsworth  (2002)  show  substantial  improvement  in 
forecast  accuracy  in  global  operational  models  over  the  past  decade,  and  present  an  extensive 
list  of  forecast  system  changes  that  may  have  contributed  to  forecast  improvement  over  the 
past  two  decades. 
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The  first  goal  of  the  present  research  is  to  compare  UGM  and  MRF  forecasts  in  more 
detail.  Our  initial  studies  of  this  problem  focused  on  a  17-case  forecast  sample  during  boreal 
winter  1993,  and  compared  research  model  forecasts  executed  at  horizontal  wavenumber  42 
triangular  truncation  with  MRF  forecasts  at  horizontal  wavenumber  62  triangular  truncation. 
Those  results  were  consistent  with  the  possibility  that  uncertainty  of  the  initial  state  repre¬ 
sented  the  greatest  limitation  to  forecast  accuracy  and  implied  that  details  of  the  model 
formulation  were  not  equally  important.  That  suggestion  is  supported  in  Fig.  1.1,  which  dis¬ 
plays  bias-corrected  anomaly  correlations  of  the  500-mb  height  forecasts  made  by  the  UGM 
and  MRF  with  analyses,  and  also  correlations  of  UGM  forecasts  with  MRF  forecasts,  averaged 
globally  (Fig  1.1a),  over  the  Southern  Hemisphere  (Figl.lb),  and  over  the  Northern  Hemi¬ 
sphere  (Figl.lc).  Bias  corrections  are  applied  by  extracting  the  average  error  of N cases  from 
the  individual  forecasts: 


N 

bias  corrected  forecast  ~  ^  -O")  ,  ( I  •  1 ) 

n=  / 

where  F*’  is  the  forecasted  variable,  O”  is  the  analyzed  variable  for  run  number  n,  and  N  is  the 
total  number  of  runs  in  the  set  of  cases  (for  Fig.  1 .1,  A''  is  17).  The  model  bias  may  also  be  cal¬ 
culated  using  a  “training  period”  applied  to  subsequent  independent  cases.  For  example,  the 
bias  can  be  computed  from  the  average  error  for  the  30  days  prior  to  a  given  date,  and  then 
applied  to  that  date.  This  “a  priori”  method  will  be  employed  in  Chapter  3  in  order  to  make 
comparison  of  forecasts  from  models  of  very  different  sophistication  more  realistic. 

The  bottom  three  curves  of  each  panel  in  Fig.  1.1,  based  on  calculations  performed  by 
MMP,  display  bias-corrected  anomaly  correlations  obtained  with  the  two  research  model  fore¬ 
casts  initialized  separately  with  National  Centers  for  Environmental  Prediction-National 
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Fig.  1.1.  Time  evolution  of  case-averaged  anomaly  correlations  of  the  500-mb  geopotential 
height,  for  the  UGM  initialized  with  NCEP-NCAR  (open  circles)  and  ECMWF  (open  squares) 
reanalyses  and  the  MRF  (triangles),  each  correlated  with  analyses,  and  for  the  UGM  initialized 
with  NCEP-NCAR  (closed  circles)  and  ECMWF  (closed  squares)  reanalyses  each  correlated 
with  the  MRF.  Correlations  are  averaged  (a)  globally,  (b)  over  the  latitude  band  from  90°  to 
20°  S,  and  (c)  over  the  latitude  band  from,  20°  to  90°  N. 
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Center  for  Atmospheric  Research  (NCEP-NCAR)  and  European  Centre  for  Medium-Range 
Weather  Forecasts  (ECMWF)  reanalyses,  and  of  the  MRF  model  initialized  with  NCEP- 
NCAR  reanalyses.  As  expected,  the  MRF  model  is  more  accurate  than  the  research  model  glo¬ 
bally  and  over  the  Northern  Hemisphere,  but  only  by  a  relatively  small  amount.  The  bias- 
corrected  UGM  lags  the  bias-corrected  MRF  model  predictions  by  only  about  12  h,  and  both 
show  similar  skill  after  seven  days  when  the  anomaly  correlations  drop  below  0.6  for  each 
model.  In  the  Southern  Hemisphere,  the  bias-corrected  UGM  initialized  with  ECMWF  data 
slightly  outperforms  the  bias-corrected  MRF  throughout  the  forecast  period.  When  calculating 
model  errors  used  for  statistics,  forecasts  from  each  model  were  not  truncated  to  the  same  res¬ 
olution.  Although  this  may  have  unfairly  penalized  the  higher-resolution  model  (MRF),  which 
has  more  finescale  structure,  the  MRF  still  outperforms  the  UGM  in  nearly  all  comparisons, 
and  without  the  bias  correction,  the  MRF  forecasts  retain  predictability  skill  24  h  longer  than 
do  the  UGM  forecasts.  The  top  two  curves  of  the  panels  in  Fig.  1.1  display  bias-corrected 
anomaly  correlations  of  two  UGM  forecasts  (initialized  with  NCEP-NCAR  and  ECMWF 
reanalyses)  validated  against  the  MRF  model  forecasts.  It  is  noteworthy  that  the  anomaly  cor¬ 
relations  of  the  UGM  forecasts  with  the  MRF  model  predictions  are  higher  than  the  anomaly 
correlation  of  either  model  with  the  analysis,  particularly  when  the  UGM  is  initialized  with 
NCEP-NCAR  reanalyses.  One  interpretation  of  these  results  is  that  the  research  model  forecast 
evolution  is  more  similar  to  that  of  the  MRF  model  than  either  model’s  forecast  evolution  is 
similar  to  the  state  of  the  atmosphere.  At  eight  days,  the  research  model  anticipates  the  MRF 
model  evolution  about  as  well  as  either  model  anticipates  the  analyzed  evolution  at  only  six 
days.  Another  interpretation  is  that  the  error  fields  of  the  MRF  and  UGM  exhibit  strong  corre¬ 
lation  through  eight  days,  implying  a  similar  error  source  for  both  models. 

The  UGM  retains  lower  horizontal  resolution  and  fewer  levels  than  the  MRF  model. 
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and  has  undergone  much  less  extensive  development  and  calibration.  The  fact  that  a  relatively 
unsophisticated  and  more  crudely  resolved  model  possesses  such  high  skill  in  anticipating  the 
behavior  of  a  much  more  sophisticated  model  implies  that  something  other  than  model  com¬ 
plexity  provides  the  fundamental  limitation  in  forecast  accuracy  for  this  set  of  17  individual 
cases  of  the  boreal  winter  of  1993.  A  possible  inference  is  that  both  models  suffer  from  the 
common  deficiencies  of  imperfect  initial  state  specification,  and  this  may  represent  the  pri¬ 
mary  forecast  limitation  for  these  cases.  This  inference  carries  important  implications  for 
deployment  of  observing  systems  and  for  forecast  model  development,  and  merits  close 
scrutiny. 

Operational  model  sophistication,  resolution,  and  initial  state  specification  have  under¬ 
gone  substantial  improvement  in  the  past  decade  (Simmons  and  Hollingsworth  2002).  The 
second  goal  of  the  current  research  is  to  evaluate  this  improvement  using  an  up-to-date  version 
of  operational  model  products  from  NCEP  and  to  compare  these  results  to  the  research  model. 
The  comparison  allows  a  crude  estimate  of  value  added  by  the  improved  observing  system  and 
data  assimilation  system  since  1993,  in  comparison  to  model  enhancements.  Chapter  2 
describes  the  datasets  and  models  used  in  the  prior  1993  comparisons,  and  more  recent  datasets 
and  models  in  real-time  comparisons.  Chapter  3  presents  updated  diagrams  similar  to  Fig.  1.1 
and  quantifies  comparisons  of  recent  research  and  operational  products.  These  results  suggest 
an  extension  of  forecast  skill  by  approximately  6  h  in  the  low-resolution  UGM  initialized  with 
boreal  winter  2003  analyses,  and  a  more  substantial  extension  of  approximately  36  h  in  a  more 
developed,  high-resolution  version  of  the  MRF  model,  due  both  to  initial  state  and  model 
improvements. 

There  is,  therefore,  evidence  for  the  important  positive  role  of  both  improved  initial 
states  as  well  as  improved  model  resolution  and  formulation.  Error  diagnostics  still  suggest  a 
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strong  similarity  of  error  structure,  consistent  with  the  possibility  that  initial  state  error  is  yet  a 
strongly  limiting  element.  Chapter  5  continues  past  studies  by  MMP  of  the  rate  of  growth  of 
the  initial  state  uncertainty.  MMP’s  experiments  were  limited  by  several  factors.  First,  they 
were  integrated  for  only  five  days,  and  it  was  not  possible  to  establish  the  timescale  on  which 
the  uncertainty  saturates.  Second,  they  did  not  demonstrate  the  relative  sensitivity  of  the  fore¬ 
casts  to  imposed  initial  uncertainty  in  comparison  to  the  sensitivity  to  model  formulation.  The 
integrations  were  also  performed  at  a  relatively  coarse  wavenumber  42  resolution,  which  could 
have  limited  error  growth. 

Therefore,  a  third  goal  of  the  present  study  is  to  perform  a  series  of  further  predictabil¬ 
ity  studies  with  the  UGM  and  compare  the  results  with  those  from  the  MRF,  with  experiments 
to  address  each  of  these  issues.  In  particular,  all  wavenumber  42  UGM  experiments  performed 
by  MMP  are  now  extended  from  five-day  to  two-week  duration;  two  sets  of  initial  condition 
sensitivity  experiments  impose  doubled  horizontal  resolution,  and  the  sensitivity  of  forecast 
evolution  to  the  model  used  is  now  compared  with  the  sensitivity  to  the  initial  state  uncertainty. 
Additionally,  UGM  cases  from  1993  and  2003  were  repeated  at  wavenumber  84  truncation  to 
determine  the  influence  of  increased  horizontal  resolution  upon  forecast  skill  and  sensitivity  to 
initial  state  uncertainty. 

A  final  goal  of  the  study  is  to  compare  results  obtained  from  the  two  primitive  equation 
models  (UGM  and  MRF)  with  those  from  a  third  model  with  a  more  general  dynamical  core. 
The  third  model  contains  an  Euler  dynamical  core  that  retains  nonhydrostatic  effects  as  well  as 
compressibility,  and  is,  therefore,  less  restricted  than  the  primitive  equation  models.  Chapter  7 
and  Appendix  A  provide  a  description  of  the  Euler  model  used  for  the  comparison,  and  indi¬ 
cate  that  there  should  be  no  systematic  differences  between  UGM  and  Euler  model  forecasts 
executed  at  present  resolutions.  Subsequent  chapters  will  show,  however,  that  inherent  model 
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differences  do  exist  and  the  choice  of  model  dynamics  may  impact  predictive  skill. 

The  study  is  organized  as  follows.  Chapter  2  presents  brief  overviews  of  the  UGM, 
MRF  model,  and  datasets.  Discussion  of  the  Euler  model  is  deferred  to  later  chapters.  Chapter 
3  describes  anomaly  correlations  of  the  UGM  predictions  in  comparison  to  predictions  by  the 
MRF  model  for  a  boreal  winter  period  in  2002/03.  Chapter  4  quantifies  the  sensitivity  of  global 
rms  estimates  of  height  and  wind  forecast  changes  produced  by  switching  between  available 
models  and  switching  between  available  initial  states,  and  compares  both  to  analyzed  and  fore¬ 
casted  changes  and  to  the  errors  of  the  changes.  Chapter  5  studies  the  time  evolution  and 
correlations  of  anomaly  correlations  during  boreal  winter  2002/03,  and  sensitivity  to  initial 
state  uncertainty  for  model  runs  at  low  resolution.  Chapter  6  illustrates  the  impact  of  doubled 
horizontal  resolution  on  predictive  skill  and  sensitivity  to  initial  state  uncertainty  for  the  UGM. 
Chapter  7  documents  the  Euler  model  used  in  the  present  study,  and  touches  on  possible  inher¬ 
ent  differences  between  that  model  and  the  UGM.  Chapter  8  describes  anomaly  correlations  of 
Euler  model  predictions,  and  compares  bias  fields  between  the  two  primitive  equation  models 
and  the  Euler  model.  Chapter  9  quantifies  the  sensitivity  of  global  rms  estimates  of  forecast 
changes  produced  by  switching  between  models  with  differing  dynamical  cores,  and  Chapter 
10  presents  overall  conclusions. 


CHAPTER  2 


DATASETS  AND  MODELS 

Datasets 

Both  NCEP,  in  collaboration  with  NCAR  (NCEP;  Kalnay  et  al.  1996),  and  the 
ECMWF  (Gibson  et  al.  1997)  have  performed  gridded  retrospective  analyses,  based  upon  all 
available  observations,  by  a  frozen  state-of-the-art  global  data  assimilation  system.  Present 
estimates  of  initial  state  uncertainty  are  obtained  from  the  difference  of  these  two  equally  cred¬ 
ible  analyses  and  are  assumed  to  be  reasonable  for  this  study.  However,  it  is  likely  that  this 
method  of  characterizing  initial  state  uncertainty  underestimates  actual  values  at  all  scales.  In 
fact,  NCEP-NCAR  reanalyses  have  been  truncated  at  wavenumber  36,  and,  to  the  degree  that 
ECMWF  and  NCEP-NCAR  use  the  same  observations,  the  difference  in  their  analyses  will 
underestimate  the  total  error.  This  limitation  to  the  study  is  further  discussed  in  Chapter  5. 

The  17  cases  previously  used  by  MMP  were  selected  for  model  initialization.  These 
start  on  1  January  1993,  and  continue  at  five-day  intervals  through  March  1993.  The  particular 
dates  are  chosen  because  eight-day  predictions  by  the  version  of  the  MRF  used  in  the  reanaly¬ 
ses  are  available  within  the  NCEP-NCAR  reanalysis  archive.  This  allows  comparison  of  the 
UGM  research  model  with  the  MRF  model,  which  is  well  documented  within  NCEP  technical 
reports. 

In  addition  to  the  earlier  1993  ensemble,  a  set  of  52  cases  was  selected  from  boreal 
winter  2002/03  from  operational  NCEP  Global  Data  Assimilation  System  (GDAS)  analyses 
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that  were  available  in  near-real  time.  These  analyses  are  archived  at  2.5°  resolution  on  26  ver¬ 
tical  levels  and  extend  from  4  January  2003  to  26  February  2003  (two  dates  were  excluded). 
The  set  of  52  cases  was  selected  based  upon  continuous  analysis  data  availability  to  compute 
bias  statistics,  allowing  an  a  priori  bias  extraction  based  upon  a  training  period  of  30-day  dura¬ 
tion  prior  to  each  UGM  and  MRP  prediction.  All  errors  and  anomaly  correlations  shown  sub¬ 
sequently  represent  bias-corrected  values  for  both  the  UGM  and  the  MRF  unless  otherwise 
stated. 

Within  the  52-case  sample  from  boreal  winter  2002/03,  a  subset  of  30  days  was  chosen 
to  initialize  the  Euler  model  runs  discussed  in  later  chapters.  Additionally,  an  even  smaller,  but 
still  representative,  subset  of  15  cases  (alternate  days  within  the  30-case  sample)  was  used  for 
high-resolution  UGM  model  initialization.  Because  the  Euler  model  time  step  is  generally  one- 
tenth  of  that  used  in  the  UGM,  and  the  high-resolution  UGM  time  step  is  one-third  that  of  the 
low-resolution  UGM,  larger  sample  sizes  were  not  possible  due  to  computational  consider¬ 
ations. 


Models 

The  UGM  is  based  upon  Galerkin  approximations  applied  separately  in  each  spatial 
dimension.  Latitude  and  vertical  structure  are  depicted  by  finite  elements  and  longitude  vari¬ 
ability  by  Fourier  series.  The  method  retains  the  high  accuracy  and  conservative  properties  of 
alternative  Galerkin  approximations,  such  as  those  used  in  the  MRF  that  are  based  upon  spher¬ 
ical  harmonic  expansions. 

The  dynamical  core  of  the  model  uses  a  hydrostatic  set  of  primitive  equations  in  which 
vorticity,  divergence,  and  thermal  fields  are  predicted  on  pressure-based  sigma  coordinates. 
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This  approach  is  similar  to  that  used  in  global  models  in  operational  centers,  with  the  exception 
of  numerical  methods  outlined  above.  Model  physical  parameterizations  of  convective  and  sta¬ 
ble  precipitation  are  similar  to  those  used  by  the  NCAR  Community  Climate  Model  1  (Bath  et 
al.  1987)  introduced  in  1987,  and  parameterizations  of  radiative  and  surface  processes  are  also 
relatively  less  advanced  and  follow  methods  that  have  been  used  by  other  models  more  than 
two  decades  ago.  Vertical  mixing  coefficients  are  calculated  from  a  low-order  turbulent  energy 
equation,  and  radiation  processes  include  cloud  radiation  interactions,  as  described  for  another 
model  by  Nicolini  et  al.  (1993).  Moist  convective  processes  use  simple  convective  adjustment, 
and  cloud  fraction  and  condensation  criteria  are  based  upon  local  relative  humidity  without 
explicit  treatment  of  cloud  microphysics.  Sea  surface  temperature  is  maintained  at  the  initial 
value,  and  land  surface  evaporation  is  extracted  from  the  daily  analysis  files.  Each  of  these  pro¬ 
cesses  is  treated  in  a  more  modem  fashion  within  the  MRF.  Present  applications  retain  20  ver¬ 
tical  levels  and  wavenumber  42  resolution.  One  set  of  experiments  was  repeated  at 
wavenumber  84  resolution. 

The  UGM  was  originally  designed  by  Paegle  (1989)  to  address  predictability  ques¬ 
tions.  It  has  been  used  to  study  impact  of  wind  data  voids  on  objective  analyses  (Paegle  and 
Horel  1991),  for  predictability  work  (Vukicevic  and  Paegle  1989;  Paegle  et  al.  1997;  Miguez- 
Macho  and  Paegle  2000,  2001),  for  Idealized  global  simulations  of  tropical-extratropical  inter¬ 
actions  (Buchmann  et  al.  1995),  to  study  orographically  forced  regional  circulations  (Nogues- 
Paegle  et  al.  1998,  Byerle  and  Paegle  2003),  and  for  initial  data  impact  investigations  of  the 
1993  “storm  of  the  centuiy”  (Miguez-Macho  and  Paegle  1999a,b). 

The  MRF  model  used  in  the  1993  intercomparisons  corresponds  to  the  version  used 
operationally  at  NCEP  until  about  1995  and  is  truncated  (triangular)  at  horizontal  wavenumber 
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62  on  27  vertical  levels.  The  acronym  for  the  MRF  changed  to  GFS  (Global  Forecasting  Sys¬ 
tem)  recently,  but  products  from  the  real-time  NCEP  forecasts  used  for  the  boreal  winter 
2002/03  cases  will  continue  to  be  referred  to  as  MRF  predictions  in  this  study.  These  forecasts 
were  executed  at  wavenumber  T254  resolution  for  the  first  3.5  days  and  subsequently  at  wave- 
number  T170  resolution  with  triangular  truncation,  on  64  vertical  levels  (S.  J.  Lord,  2003,  per¬ 
sonal  communication).  The  MRF  is  a  spectral  model,  based  upon  a  set  of  primitive  equations 
similar  to  those  used  in  the  UGM. 

The  version  of  the  Euler  model  used  in  the  present  study  is  relatively  new,  and  has  not 
been  previously  documented.  Therefore,  detailed  discussion  of  this  model  will  be  deferred  to  a 
separate  chapter  and  appendix. 


CHAPTER  3 


ANOMALY  CORRELATIONS 

Figure  3.1  compares  500-mb  geopotential  height  anomaly  correlations  for  the  recent 
52-case  real-time  forecasts  of  boreal  winter  2002/03  with  those  for  the  1993  experiments,  aver¬ 
aged  globally  (Fig.  3.1a)  and  over  the  boreal  winter  (northern)  hemisphere  (Fig.  3.1b).  The 
anomaly  correlations  of  forecasts  calibrated  against  observations  (in  this  case,  analyses)  are 
computed  as 


anom  corr  =  , 


(3.1) 


where  S^i  represents  the  climatological  value  of  S  and  S^i,  is  the  analyzed  value,  in  this  case 
taken  from  GDAS  analyses.  Here  ( )  indicates  a  global  or  hemispheric  average,  and  the  50-yr 
(1951-2000)  monthly  average  derived  from  the  NCEP-NCAR  reanalyses  is  defined  as  the  cli¬ 
matology  for  this  purpose;  S  represents  the  forecasted  variable  from  the  UGM  and  MRF.  An  a 
priori  bias  correction  is  applied  to  forecast  variables  using  biases  calculated  from  the  30  daily 
forecasts  available  immediately  prior  to  the  day  of  interest.  The  anomaly  correlations  of  UGM 
forecasts  calibrated  against  MRF  forecasts  are  given  by  a  similar  formula,  with  S  replaced  by 
the  UGM  prediction  and  5^^  by  the  MRF  prediction.  Anomaly  correlations  of  the  UGM  fore¬ 
casts  with  the  MRF  model  predictions  exceed  anomaly  correlations  of  UGM  forecasts  with 
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Fig.  3.1.  Time  evolution  of  case-averaged  anomaly  correlations  of  the  500-mb  geopotential 
height,  for  the  1993  cases  of  the  UGM  (open  circles)  and  the  MRF  (open  squares),  and  for  the 
boreal  winter  2002/03  cases  of  the  UGM  (closed  circles)  and  MRF  (closed  squares),  each  cor¬ 
related  with  analyses;  and  for  the  boreal  winter  2002/03  cases  of  the  UGM  correlated  with  the 
MRF  (closed  triangles.  Curves  are  averaged  (a)  globally  and  (b)  over  the  latitude  band  from 
20°  to90°N. 


Fig.  3.1.  (continued) 
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analyses  when  errors  of  the  UGM  possess  similar  structure  to  errors  of  the  MRF. 

The  curves  with  solid  squares  in  Fig.  3.1  represent  anomaly  correlations  of  MRF  pre¬ 
dictions  with  analyses  for  the  52  selected  cases  from  boreal  winter  2002/03.  Substantial 
improvement  is  seen  with  respect  to  similar  results  for  the  1993  ensemble  (curves  with  open 
squares).  Bias-corrected  MRF  anomaly  correlations  for  boreal  winter  2002/03  drop  below  0.6 
after  approximately  192  h,  compared  to  approximately  156  h  for  the  1993  cases,  corresponding 
to  a  36  h  gain  in  predictive  skill  at  this  level  by  the  MRF. 

The  UGM  also  produces  more  skillful  forecasts  in  2003  than  it  does  in  1993.  In  the 
boreal  winter  hemisphere  (Fig.  3.1b),  UGM  forecasts  produce  anomaly  correlations  of  0.6  out 
to  approximately  162  h  (curve  with  solid  circles)  and  represent  6  h  longer  skill  than  do  those 
for  the  1 993  ensemble  (open  circles). 

The  MRF  has  clearly  benefited  more  from  advances  made  at  NCEP  during  the  past 
decade  than  has  the  UGM.  In  the  case  of  the  MRF,  these  advances  include  substantially  higher 
horizontal  and  vertical  resolution  in  addition  to  improved  initial  state  analyses.  The  only 
improvement  in  the  UGM  is  the  availability  of  more  modern  analyses  of  the  initial  state,  since 
UGM  resolution  and  physical  formulation  are  frozen  between  the  1993  and  more  recent  appli¬ 
cations  displayed  in  Fig.  3.1. 

A  superficial  examination  of  Fig.  3.1  suggests  that  observation  and/or  analysis 
improvements  have  produced  only  a  6  h  improvement  in  forecasts  over  the  past  decade,  since 
the  frozen  UGM  exhibits  only  a  6  h  forecast  improvement  in  that  period.  It  is  unclear,  however, 
how  much  of  the  additional  30  h  of  forecast  improvement  in  the  MRF  is  due  to  initial  state 
advances  relative  to  model  advances  over  the  past  decade,  or  a  relationship  between  the  two, 
since  the  first  guess  is  such  an  important  component  of  the  data  assimilation  cycle.  A  substan- 
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tial  portion  of  that  model’s  improvement  may  be  due  to  the  improved  specification  of  the  initial 
state  in  addition  to  model  advances.  It  is  also  possible  that  boreal  winter  2003  was  more  pre¬ 
dictable  than  boreal  winter  1993  and  this  may  improve  recent  skill  scores  in  all  models. 

The  curves  with  triangles  in  Fig.  3.1  represent  anomaly  correlations  of  the  UGM  fore¬ 
casts  with  forecasts  made  by  the  MRF  for  the  boreal  winter  2002/03  cases.  This  curve  is  only 
slightly  below  the  anomaly  correlation  of  the  MRF  forecasts  with  respect  to  analyses  (curve 
with  solid  squares)  in  the  boreal  winter  hemisphere  (Fig.  3.1b).  Consequently,  the  low-resolu¬ 
tion  UGM  forecasts  have  almost  as  much  skill  in  anticipating  the  MRF  forecasts  as  the  latter 
have  in  anticipating  the  state  of  the  atmosphere.  This  implies  a  correlation  of  error  patterns  of 
forecasts  made  by  the  UGM  and  by  the  operational  model.  The  similarity  of  the  error  patterns 
of  two  very  different  models  suggests  common  limitations,  including  limitations  associated 
with  use  of  the  same  initial  state. 

Figure  3.2  shows  an  example  of  the  bias-corrected  error  fields  for  the  500-mb  meridi¬ 
onal  wind  for  day  14  (336  h)  of  forecasts  made  by  the  UGM  (Fig.  3.2a)  and  MRF  (Fig.  3.2b). 
Accurate  prediction  of  the  meridional  flow  is  important  for  accurate  depiction  of  wave  ampli¬ 
tude  and  phase,  and  this  variable  is  more  representative  of  synoptic-scale  features  than  geopo¬ 
tential  height,  which  carries  the  most  amplitude  on  global  scales.  Error  field  structural 
similarities  in  the  500-mb  meridional  flow  are  particularly  evident  in  northern  Canada,  over 
the  Eurasian  subcontinent,  and  the  eastern  coast  of  Asia  in  the  Northern  Hemisphere.  In  the 
Southern  Hemisphere,  two  strong  error  field  couplets  exist  off  the  southwestern  coasts  of 
Africa  and  Australia  in  both  forecasts.  Correlation  coefficients  (discussed  below)  for  the  bias- 
corrected  error  fields  shown  in  Fig.  3.2  are  0.61  for  global  and  Northern  Hemisphere  averages, 
and  0.62  for  the  Southern  Hemisphere  average.  Similar  agreement  in  error  patterns  of  the  two 


Fig.  3.2.  Bias-corrected  error  fields  for  the  500-mb  meridional  wind  forecast  at  hour  336,  ini¬ 
tialized  at  0000  UTC  on  26  Feb  2003,  for  the  (a)  UGM  and  (b)  MRF  model.  Contour  interval  is 
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model  forecasts  is  evident  on  many  other  days. 

These  error  fields  are  computed  after  removing  from  each  model  the  bias  calculated 
from  the  previous  30  days,  shown  in  Fig.  3.3.  The  bias  fields  for  the  500-mb  meridional  wind 
also  display  remarkable  similarity  between  the  UGM  (Fig.  3.3a)  and  MRF  (Fig.  3.3b)  at  day  14 
(336  h),  despite  the  difference  in  model  complexity  and  resolution,  implying  some  common 
deficiency.  Correlation  coefficients  for  the  30-day  bias  fields  shown  in  Fig.  3.3  are  0.74,  0.72, 
and  0.81  for  the  global  average.  Northern  Hemisphere,  and  Southern  Hemisphere,  respectively. 
These  clearly  exceed  the  asymptotic  correlation  limit  of  0.5,  discussed  in  appendix  B.  Further¬ 
more,  the  overall  magnitudes  of  the  30-day-average  meridional  flow  biases  from  the  models 
(Figs.  3.3)  are  comparable  to  the  magnitude  of  the  30-day-average  500-mb  analysis  (Fig.  3.4) 
for  meridional  wind  (Fig.  3.4a).  The  model  biases,  therefore,  are  not  negligible. 

A  comparison  of  the  departures  from  the  zonal  mean  for  the  30-day  average  500-mb 
heights  from  the  reanalysis  (Fig.  3.4b)  with  the  departures  from  the  zonal  mean  of  30-day  aver¬ 
ages  of  500-mb  forecasted  heights  (Fig.  3.5)  from  the  UGM  (Fig.  3.5a)  and  MRF  (Fig.  3.5b) 
also  shows  a  stronger  similarity  between  model  forecasts  than  between  either  model  and  the 
GDAS  analysis.  Both  the  UGM  (Fig.  3.5a)  and  MRF  (Fig.  3.5b)  tend  to  forecast  a  ridge  over 
the  northern  Rockies,  while  the  reanalysis  (Fig.  3.4b)  indicates  that  the  actual  ridge  placement 
is  off  the  western  coast  of  the  United  States  and  Canada.  Additionally,  neither  model  forecasts 
the  blocking  pattern  over  Europe,  and  both  models  underforecast  the  amplitude  of  most  of  the 
wave  structure  in  the  500-mb  pattern. 

To  further  quantify  the  similarities.  Figs.  3.6  and  3.7  display  correlation  coefficients  of 
the  error  and  bias  of  the  500-mb  meridional  flow  forecast  by  the  UGM  correlated  with  the  error 
and  bias  of  the  same  variable  forecast  by  the  MRF: 
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Fig.  3.6.  Time  evolution  of  case-averaged  correlation  coefficients  (not  bias  corrected)  of  the 
500-mb  meridional  wind  error  for  30  days  preceding  26  Feb  2003,  averaged  globally  (solid), 
averaged  over  the  latitude  band  90°  -20°  S  (open  circles),  and  averaged  over  the  latitude  band 
20°-90°N  (open  squares). 
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Fig.  3.7.  Time  evolution  of  correlation  coefficients  of  the  500-mb  meridional  wind  bias,  com¬ 
puted  over  the  30  days  preceding  26  Feb  2003,  averaged  globally  (solid),  averaged  over  the 
latitude  band  90°  -20°  S  (open  circles),  and  averaged  over  the  latitude  band  20°  -90  °N  (open 


corr  coeff = 


(3.2) 


where  5^/(5;^  represents  either  the  error  or  bias  of  the  500-nib  flow  forecasted  by  the  UGM,  and 
Smrf  is  either  the  error  or  bias  of  the  500-mb  flow  forecasted  by  the  MRF.  Again,  Q  indicates 
a  global  or  hemispheric  average.  This  measure  of  correlation  asymptotes  to  a  nonzero  value  at 
long  times  (see  appendix  B),  when  predictability  is  lost  due  to  initial  state  uncertainty  and  the 
chaotic  nature  of  the  atmosphere.  The  time  evolution  of  the  case-averaged  error  field  correla¬ 
tion  coefficients  for  a  30-case  sample  in  boreal  winter  2003  is  shown  in  Fig.  3.6.  The  fact  that 
the  correlations  of  errors  continue  to  rise  throughout  the  14-day  forecasts  despite  initial  data 
uncertainties  and  chaos  may  be  explained  by  the  growing  similarities  in  500-mb  meridional 
wind  model  bias  between  the  UGM  and  MRF,  shown  by  the  correlation  coefficients  in  Fig.  3.7. 
As  yet,  however,  it  is  unclear  to  what  the  similarity  in  the  bias  structure  can  be  attributed. 

Given  the  similarity  of  errors  in  models  of  highly  differing  complexity,  it  is  tempting  to 
speculate  that  a  fundamental  limitation  of  forecast  accuracy  is  the  uncertainty  of  the  initial 
state,  a  problem  common  to  all  global  models  initialized  with  the  same  or  similar  data,  or 
inherent  errors  in  model  physics  or  parameterizations  common  to  the  global  models  in  ques¬ 
tion.  The  next  section  describes  other  measures  of  the  forecast  sensitivity  to  model  and  initial 
state  changes. 


CHAPTER  4 


GLOBAL  ROOT-MEAN-SQUARE  FORECAST 
AND  ANALYZED  EVOLUTION 

Although  UGM  and  MRF  model  error  fields  display  substantial  correlation,  the  fore¬ 
casts  also  exhibit  differences  and  these  are  most  easily  summarized  in  terms  of  their  global  rms 
values: 


- j  1/2 

rms  =  ((P-a))  ,  (4.1) 

where  p  and  a  are  some  time-specific  variable  (as  described  below)  and  (  )  indicates  a  global 
average.  Figure  4.1  displays  global  rms  changes  of  the  analyzed  height/wind  (P  is  the  ana¬ 
lyzed  height/wind  at  time  t;  a  is  the  analyzed  height/wind  at  initial  time)  as  well  as  the  global 
rms  magnitude  of  the  forecast  change  in  height/wind  (P  is  the  forecast  value  at  time  t;  a  is  the 
analyzed  value  at  initial  time).  It  also  displays  the  model  forecast  error  for  height/wind  (P  is 
the  forecast  value  at  time  t;  a  is  the  analyzed  value  at  time  t),  and  the  sensitivity  of  the 
height/wind  forecast  to  the  model  used  (p  is  the  UGM  forecast  value  at  time  t;  a  is  the  MRF 
forecast  value  at  time  t).  Finally,  the  figure  displays  the  sensitivity  of  height/wind  forecasts  to 
model  initial  condition  changes  (P  is  the  NCEP-NCAR-initialized  UGM  forecast  value  at  time 
t;  a  is  the  ECMWF-initialized  forecast  value  at  time  t). 

The  top  (solid)  curves  of  Fig.  4.1  depict  the  global  rms  evolution  of  the  1993  analyzed 


Fig.  4.1.  Time  evolution  of  selected  global,  case-averaged  rms  variables,  for  the  1993  cases. 
The  top  curve  (solid)  depicts  the  differences  of  the  analyzed  variable  at  the  time  indicated 
along  the  abscissa  minus  the  value  at  0  h.  The  second  curve  (open  circles)  depicts  the  differ¬ 
ences  of  the  MRF  forecast  variable  at  the  indicated  time  minus  the  analysis  at  0  h.  The  third 
curve  (closed  circles)  depicts  the  differences  of  the  UGM  (initialized  with  NCEP-NCAR)  vari¬ 
able  at  the  indicated  time  minus  the  analysis  at  0  h.  The  fourth  curve  (open  squares)  depicts  the 
rms  error  of  the  variable  forecast  by  the  MRP.  The  fifth  curve  (closed  squares)  depicts  the  rms 
error  of  the  variable  forecast  by  the  UGM  (initialized  with  NCEP-NCAR).  The  sixth  curve 
(xs)  depicts  the  difference  of  the  UGM  (initialized  with  NCEP-NCAR)  forecast  variable  minus 
the  MRF  forecast  variable,  and  measures  the  sensitivity  of  the  forecasts  to  the  model  used.  The 
seventh  and  bottom  curve  (diamonds)  depicts  the  difference  of  the  UGM  forecast  variable  ini¬ 
tialized  with  NCEP-NCAR  reanalyses  minus  that  of  the  UGM  initialized  with  ECMWF 
reanalyses,  and  measures  the  sensitivity  of  the  forecasts  to  the  initial  dataset  used.  Selected 
variables  are  500-mb  (a)  height  and  (b)  meridional  wind. 


Fig.  4.1.  (continued) 
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differences  of  the  500-mb  height  (Fig.  4.1a)  and  meridional  wind  (Fig.  4.1b)  fields,  respec¬ 
tively,  at  the  time  indicated  along  the  abscissa  minus  the  initial  value.  This  curve  is  an  approx¬ 
imate  measure  of  the  observed  change  in  the  "signal"  associated  with  each  field  (i.e.,  the  error 
of  a  persistence  forecast)  and  represents  the  average  for  the  17  cases.  The  curve  in  each  dia¬ 
gram  immediately  below  this  (open  circles),  and  stopping  at  192  h,  displays  the  rms  evolution 
of  the  respective  field  from  initial  to  forecast  time  for  the  eight-day  bias-corrected  MRF  fore¬ 
casts.  The  fact  that  these  curves  are  slightly  lower  than  the  topmost  curves  suggests  that  the 
MRF  model  maintains  slightly  less  variability  than  is  analyzed  in  the  atmosphere. 

This  deficiency  is  even  more  pronounced  in  the  low-resolution  UGM,  for  which  the 
bias-corrected  forecast  signal  is  displayed  in  the  third  curve  from  the  top  (closed  circles)  in 
Figs.  4.1a  and  4.1b.  The  forecasted  signal  of  the  UGM  is  approximately  10%-20%  weaker  than 
the  analyzed  signal.  The  problem  is  particularly  evident  in  the  geopotential  height  evolution, 
for  which  the  rms  change  from  the  initial  state  continues  to  grow  past  day  14  (336  h),  while 
this  measure  of  evolution  grows  very  little  past  day  seven  (1 68  h)  in  the  UGM  (even  though  the 
UGM  predicted  fields  continue  to  change  after  day  seven). 

The  fourth  curve  from  the  top  (open  squares),  stopping  at  192  h  depicts  the  globally 
computed  rms  error  of  the  500-mb  geopotential  height  (Fig.  4.1a)  and  meridional  wind  (Fig. 
4.1b),  for  the  MRF.  The  fifth  curve  from  the  top  (closed  squares)  depicts  the  globally  computed 
rms  error  of  the  500-mb  geopotential  height  (Fig.  4.1a)  and  meridional  wind  (Fig.  4.1b),  for  the 
NCEP-NCAR  initialized  UGM.  These  curves  asymptote  toward  the  forecasted  signal  after  day 
8  and  the  relative  error  growth  is  slowest  for  the  height  field  and  faster  for  the  meridional  wind 
field,  as  expected. 

The  bottom  two  curves  in  Figs.  4.1a  and  4.1b  display  the  sensitivity  of  the  forecast  to 
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the  model  used  (second  from  the  bottom,  xs)  and  to  the  initial  dataset  used  (diamonds).  The 
curves  depicting  the  sensitivity  to  the  model  used  are  obtained  by  differencing  the  forecast  pro¬ 
duced  by  the  MRF  model  from  that  produced  by  the  low-resolution  UGM  and  computing  the 
globally  averaged  rms  of  the  result  for  the  17  cases.  Through  the  eight-day  period  for  which 
this  difference  field  can  be  calculated,  this  measure  of  sensitivity  to  the  model  is  less  than  the 
forecast  error  of  either  the  UGM  or  the  MRF  predictions  in  each  of  the  displayed  fields.  This 
supports  the  conjecture  that  forecasts  by  different  models  begun  with  the  same  initial  condi¬ 
tions  sometimes  resemble  each  other  more  closely  than  they  resemble  reality. 

The  bottom  curves  of  Figs.  4.1a  and  4.1b  show  sensitivity  to  a  data  switch  from  NCEP- 
NCAR  to  ECMWF  reanalyses  for  initial  data.  The  fact  that  this  curve  lies  well  below  the  two 
curves  directly  above  it  (forecast  error  and  model  sensitivity)  may  imply  that  initial  data  uncer¬ 
tainty  is  not  a  very  important  problem  in  the  1 993  cases.  This  inference  does  not  have  strong 
support,  since  certain  other  features  of  the  forecast  evolution  imply  sensitivity  to  the  initial 
data,  as  discussed  in  the  previous  section.  It  is  also  likely  that  the  present  estimate  of  initial 
uncertainty  (from  the  difference  of  ECMWF  and  NCEP-NCAR  analyses  of  essentially  the 
same  observations)  underestimates  the  actual  initial  state  errors.  Additionally,  the  rate  of 
increase  of  this  curve  is  related  to  model  resolution,  and  may  be  a  function  of  model  complex¬ 
ity  as  well.  These  aspects  are  discussed  in  more  detail  in  Chapter  5. 

Figures  4.2  displays  the  same  data  for  the  52  selected  cases  of  boreal  winter  2002/03, 
with  all  forecasts  initialized  with  GDAS  analyses.  In  general,  most  of  the  discussion  for  the 
1993  cases  also  applies  here,  although  the  model  forecast  error  curves  have  smaller  values  in 
2003  than  in  1993  for  both  the  UGM  (closed  squares)  and  the  MRF  (open  squares),  while  the 
sensitivity  to  the  model  (xs)  has  increased  slightly.  For  boreal  winter  2002/03,  the  sensitivity  of 
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Fig.  4.2.  Time  evolution  of  selected  global,  case-averaged  rms  variables,  for  the  2003  cases. 
The  top  curve  (solid)  depicts  the  differences  of  the  analyzed  variable  at  the  time  indicated 
along  the  abscissa  minus  the  value  at  0  h.  The  second  curve  (open  circles)  depicts  the  differ¬ 
ences  of  the  MRP  forecast  variable  at  the  indicated  time  minus  the  analysis  at  0  h.  The  third 
curve  (closed  circles)  depicts  the  differences  of  the  UGM  variable  at  the  indicated  time  minus 
the  analysis  at  0  h.  The  fourth  curve  (open  squares)  depicts  the  rms  error  of  the  variable  fore¬ 
cast  by  the  MRP.  The  fifth  curve  (closed  squares)  depicts  the  rms  error  of  the  variable  forecast 
by  the  UGM.  The  sixth  curve  (xs)  depicts  the  difference  of  the  UGM  forecast  variable  minus 
the  MRP  forecast  variable,  and  measures  the  sensitivity  of  the  forecasts  to  the  model  used. 
Selected  variables  are  500-mb  (a)  height  and  (b)  meridional  wind. 
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the  forecast  to  the  selected  model  is  about  as  great  as  the  error  of  the  MRF  model,  and  only 
slightly  smaller  than  the  error  of  the  UGM  through  the  first  week  of  prediction.  This  result  and 
the  substantial  positive  correlation  between  UGM-  and  MRP-predicted  error  fields  discussed 
in  the  last  section  suggest  that  the  two  models  produce  errors  that  have  similar  patterns  that  dif¬ 
fer  in  magnitude.  This  “distinct  tendency  for  errors  in  each  model  to  be  highly  correlated  with 
one  another,”  though  remarkable,  is  not  new  and  was  noted  in  forecast  intercomparisons  over 
two  decades  ago  (Baumhefner  and  Downey  1978). 


CHAPTER  5 


INITIAL  STATE  ERROR  EVOLUTION 

The  similarity  of  forecast  errors  between  the  UGM  and  GFS  is  reflected  in  positive  cor¬ 
relations  in  the  skill  scores  of  the  two  models.  Figure  5.1  displays  a  52-day  time  series  from 
boreal  winter  2002/03  of  globally-averaged  500-mb  height  anomaly  correlations  of  UGM 
(curves  with  circles)  and  MRF  (curves  with  squares)  forecasts,  validated  against  NCEP’s 
GDAS  analyses,  and  the  anomaly  correlation  of  UGM  forecasts  validated  against  MRF  fore¬ 
casts  (triangles).  Results  are  shown  separately  for  forecast  days  1,  3,  5,  7,  9  and  1 1  (Figs.  5.1a, 
b,  c,  d,  e,  and  f,  respectively). 

There  is  a  tendency  for  the  anomaly  correlations  of  the  UGM  to  show  similar  time  vari¬ 
ation  to  those  of  the  MRF.  Skill  scores  of  the  two  models  frequently  increase  and  decrease  in 
tandem,  resulting  in  higher  values  of  model-to-model  correlations  during  these  time  periods. 
This  trend  is  quantified  in  Fig.  5.2,  which  depicts  the  variation  with  forecast  day  of  the  correla¬ 
tion  coefficient  between  the  average  of  each  of  the  two  solid  curves  in  Fig.  5.1.  The  correlation 
coefficient  peaks  on  day  one  of  the  forecasts,  reduces  to  a  minimum  after  about  four  days,  and 
then  increases  in  the  second  week  of  the  forecasts.  Although  actual  correlation  values  are  low, 
the  qualitative  behavior  is  consistent  with  the  likelihood  that  both  models  are  similarly  limited 
by  initial  condition  uncertainties  early  in  their  forecasts,  and  by  similar  bias  evolutions  at 
longer  times. 

Figure  5.1  suggests  substantial  episodic  variability  in  predictive  skill  after  one  week  of 
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Fig.  5.1.  Time  series  of  globally-averaged  anomaly  correlations  of  the  500-mb  geopotential 
height,  for  the  UGM  initialized  with  GDAS  analyses  (circles),  and  the  MRF  (squares),  each 
correlated  with  analyses;  and  for  the  UGM  correlated  with  the  MRF  (triangles).  Correlations 
are  shown  at  (a)  forecast  hour  24,  (b)  forecast  hour  72,  (c)  forecast  hour  120,  (d)  forecast  hour 
168,  (e)  forecast  hour  216,  and  (f)  forecast  hour  264. 
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Fig.  5.2.  Time  evolution  of  correlation  coefficient  of  boreal  winter  2002/03  52-case  average 
bias  corrected  500-mb  height  UGM  anomaly  correlation  to  MRF  anomaly  correlation.  Time 
scale  on  the  x-axis  is  forecast  day. 
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prediction.  The  MRF  exhibits  skill  scores  around  0.8  on  several  seven-day  forecasts  (Fig. 
5. Id),  and  skill  scores  below  0.6  for  seven-day  forecasts  during  other  periods.  It  would  be  use¬ 
ful  to  determine  whether  such  highly  variable  behavior  is  due  mostly  to  fluctuations  in  the 
accuracy  of  initial  state  specification,  or  due  to  the  possibility  that  different  atmospheric  con¬ 
figurations  are  more  difficult  to  predict  than  others.  This  important  area  of  research  lies  outside 
the  scope  of  the  present  study. 

The  curves  with  triangles  in  Fig  5.1  also  reflect  similarity  of  forecast  errors  between  the 
two  models.  These  curves  represent  the  similarity  of  the  UGM  predictions  to  the  evolution  of 
the  MRF.  When  the  values  of  the  model-to-model  correlations  (triangles)  exceed  those  of  the 
UGM  (circles),  the  UGM  and  MRF  forecasts  are  more  similar  to  each  other  than  the  UGM 
forecast  is  similar  to  the  atmosphere.  This  is  almost  always  the  case,  as  shown  in  Fig  5.1.  On 
those  occasions  that  the  model-to-model  correlations  (triangles)  lies  above  the  MRF  correla¬ 
tions  (squares),  the  MRF  predictions  are  more  similar  to  the  UGM  predictions  than  to  the 
atmosphere.  This  occurs  23  different  times  during  the  52  forecasts  available  for  day  11.  These 
results  suggest  that  the  errors  of  the  two  models  are  similar,  not  only  for  the  average  of  many 
forecasts  as  shown  in  Chapter  2,  but  also  for  many  individual  forecasts. 

Once  again,  results  are  consistent  with  a  common  error  source.  One  shared  source  of 
errors  is  the  initial  state,  which  is  the  same  in  all  of  the  predictions  reflected  in  Figs.  5.1  and 
5.2.  The  question  of  which  part  of  the  initial  state  uncertainty  may  provoke  the  greatest  fore¬ 
cast  sensitivity  is  now  discussed. 

MMP  demonstrated  that  the  influence  of  initial  state  uncertainty  is  far  from  saturation 
after  five  days  of  simulation  by  the  Utah  model.  They  studied  forecast  sensitivity  to  initial  state 
uncertainty  due  to  different  complementary  wave  groups  in  the  total  initial  uncertainty.  The 
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uncertainty  was  estimated  from  the  difference  of  two  equally  credible  analyses  of  the  atmo¬ 
sphere,  provided  by  NCEP-NCAR  (T36  truncation)  and  ECMWF  reanalyses  of  the  same  state 
during  boreal  winter  1993.  These  analyses,  which  use  the  same  observations,  are  very  similar 
over  land  areas  of  good  observational  coverage  (see  Fig.  5.3  for  500-mb  height  and  wind  anal¬ 
ysis  differences),  but  they  differ  more  substantially  in  regions  of  poor  observations,  such  as  the 
southern  oceans. 

In  qualitative  terms,  the  spatial  structures  of  analysis  differences  displayed  in  Fig.  5.3, 
represent  the  expected  geographical  distribution  of  the  initial  state  uncertainty.  In  particular, 
the  analysis  differences  are  small  over  well-observed  continents  and  larger  over  poorly 
observed  regions  of  the  globe.  For  this  study,  it  is  assumed  this  difference  field  provides  a  rea¬ 
sonable  estimate  of  the  actual  observation  uncertainty.  The  influence  of  this  level  of  initial 
state  uncertainty  is  determined  by  repeating  separate  forecasts  made  by  the  low-resolution 
UGM  initialized  with  NCEP-NCAR  and  with  ECMWF  reanalyses,  and  studying  the  evolving 
difference  fields  of  the  predictions.  Following  MMP,  sequences  of  experimental  cases  are  per¬ 
formed  in  which  separate  spectral  wave  groups  are  modified  in  the  initial  data,  and  results  are 
normalized  by  dividing  by  the  variance  of  the  predicted  difference  produced  by  initial  state 
modifications  of  the  entire  spectrum.  Present  results  emphasize  the  response  in  the  meridional 
wind  since  this  effectively  represents  the  positioning  and  magnitude  of  synoptic-scale  waves. 
This  quantity  also  better  represents  the  variability  in  the  Tropics  than  does  the  height  field. 

The  results  are  displayed  in  Fig.  5.4,  which  depicts  normalized  forecast  responses  to 
uncertainties  retained  in  different  portions  of  the  global  spectrum  for  the  1993  cases.  If  the 
forecasts  were  extended  over  a  sufficiently  long  period,  error  growth  would  saturate,  and  all 
the  curves  of  Figs.  5.4a  and  5.4b  would  asymptote  to  1.0,  and  the  curves  of  Fig.  5.4c  would 
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(a)  time  (h) 


Fig.  5.4.  Time  evolution  of  global  variances  at  a  =  0.53,  of  forecast  meridional  flow  differ¬ 
ences  produced  by  initial  data  changes  from  NCEP-NCAR  to  ECMWF  reanalyses  for  the  17 
cases  from  boreal  winter  1993.  (a)  Curves  show  the  impact  of  changing  the  initial  data  only  for 
global  wavenumbers  0-5  (open  circles),  0-10  (closed  circles),  0-15  (open  squares),  and  0-20 
(closed  squares),  (b)  Curves  display  the  impact  of  changing  the  initial  data  only  for  the  com¬ 
plementary  wavenumbers  6-42  (open  circles),  11-42  (closed  circles),  16-42  (open  squares), 
and  21-42  (closed  squares).  Results  in  both  (a)  and  (b)  have  been  normalized  by  dividing  by 
the  variance  of  the  forecast  meridional  flow  differences  produced  by  the  initial  state  change  of 
the  entire  spectrum.  Sums  of  complementary  wave  groups  in  (a)  and  (b)  are  displayed  in  (c). 
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Fig.  5.4.  (continued) 


normalized  variance 
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Fig.  5.4.  (continued) 
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then  asymptote  to  2.0. 

MMP  produced  the  results  of  Figs.  5.4a  and  5.4b  out  to  five  days  (120  h),  and  noted 
that  sums  of  normalized  uncertainty  growth  for  complementary  wave  groups  approximately 
equal  1.0  rather  than  2.0.  They  concluded  that  the  first  five  days  of  prediction  by  the  Utah 
model  were  characterized  by  slow,  linear  error  growth,  which  is  far  from  saturation  through 
five  days.  Within  this  regime  of  unsaturated  error  growth,  the  effect  of  changing  long-wave 
components  of  the  uncertainty  spectrum  produces  an  error  contribution  that  grows  more  rap¬ 
idly  in  a  relative  sense  than  does  the  effect  of  changing  shorter  waves  of  the  initial  state.  This 
can  be  deduced  from  the  initial  upward  trends  of  the  curves  in  Fig.  5.4a  and  the  initial  down¬ 
ward  trends  in  the  curves  in  Fig.  5.4b. 

The  present  extension  to  14-day  predictions  provides  important  modifications  of  these 
conclusions.  In  particular,  pronounced  growth  of  the  sum  of  the  normalized  uncertainty  contri¬ 
butions  in  Fig.  5.4c  occurs  after  120  h,  when  the  sums  are  on  the  order  of  1.1  or  less.  By  336  h, 
the  sums  exhibited  in  Fig.  5.4c  range  from  1.5  to  nearly  1.7  suggesting  that  initial  state  errors 
have  effectively  spread  across  the  full  spectrum  resolved  in  the  present  integrations.  Even  for 
this  extended  prediction  the  longer-wave  uncertainty  explains  at  least  as  much,  or  more,  of  the 
total  error  growth  than  does  uncertainty  in  the  shorter  waves.  Waves  0-15,  for  example, 
account  for  less  than  40%  of  the  initial  uncertainty  in  the  present  experiments  in  Fig.  5.4a,  but 
produce  about  50%  more  sensitivity  at  14  days  (336  h)  than  do  waves  16-42  (Fig.  5.4b). 

The  sums  of  the  relative  uncertainty  in  different  spectral  groups  increase  toward 
approximately  1.5-1. 7  in  Fig.  5.4c.  These  curves  should  asymptote  to  2.0  in  the  limit  when 
errors  associated  with  the  chaotic  nature  of  the  atmosphere  would  finally  saturate.  Simmons 
and  Hollingsworth  (2002)  suggest  that  the  ECMWF  model  is  close  to  error  saturation  after 
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approximately  21  days.  They  base  this  conclusion  upon  the  rate  at  which  model  solutions  ini¬ 
tialized  on  consecutive  days  diverge.  The  present  experiments  are  integrated  to  only  14  days 
and  use  a  different  method  to  study  error  saturation,  but  an  extrapolation  of  the  curves  in  Fig. 
5.4c  from  the  second  week  shown  toward  21  days  appears  to  be  consistent  with  Simmons  and 
Hollingsworth’s  (2002)  conclusion. 

The  most  important  limitations  to  the  present  conclusions  are  that  they  are  made  using 
a  forecast  model  truncated  at  relatively  low  resolution,  and  that  the  differences  of  NCEP- 
NCAR  and  ECMWF  analyses,  which  use  the  same  observations,  may  not  adequately  reflect 
the  actual  uncertainty  of  the  initial  state.  Higher-resolution  models  generally  allow  larger 
uncertainty  growth  because  model  diffusion  coefficients  need  to  increase  with  coarser  resolu¬ 
tion  to  control  spectral  blocking  and  to  limit  the  accumulation  of  energy  at  the  shortest  resolv¬ 
able  scales  due  to  nonlinear  energy  cascade. 

It  is  likely  that  the  present  method  of  initial  state  uncertainty  specification  underesti¬ 
mates  the  actual  uncertainty  in  both  large  and  small  scales.  In  particular,  neither  analysis  con¬ 
tains  much  amplitude  in  higher-wavenumber,  smaller  scale  components  (wavelengths  on  the 
order  of  1000  km)  and  their  differences  may  systematically  underestimate  the  impact  of  errors 
at  these  scales. 

White  et  al.  (1999)  estimated  “errors”  of  the  MRF  initial  state  over  the  intermountain 
west  of  the  United  States  in  a  region  of  good  radiosonde  density.  Table  2  of  that  study  suggests 
that  initial  rms  error  values  for  the  wind  exceed  4  m  s'*  in  the  mid  troposphere  in  a  region  of 
reliable  observations.  These  “uncertainties”  are  estimated  from  the  fit  of  the  initial  analyses 
interpolated  to  radiosonde  observation  sites  using  a  horizontal  grid  corresponding  to  approxi¬ 
mately  global  wavenumber  100.  Presumably,  the  uncertainties  are  even  larger  in  regions  of 
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sparse  observations.  The  initial  data  uncertainties  (Fig.  5.3)  used  in  the  present  experiments  are 
smaller  than  4  m  s"'  over  most  of  the  globe,  and  impose  smaller  initial  state  uncertainties  than 
were  found  by  White  et  al.  (1999)  in  data-rich  regions. 

In  order  to  address  these  issues,  three  other  sets  of  experiments  were  conducted.  In  one, 
the  initial  state  uncertainty  estimated  from  the  ECMWF  -  NCEP-NCAR  forecast  differences 
was  retained  only  in  wave  groups  16-42,  as  in  the  earlier  experiment  sequence,  but  the  ampli¬ 
tude  of  this  difference  field  was  arbitrarily  doubled.  A  similar  approach  was  taken  by  doubling 
the  initial  state  uncertainty  only  in  waves  0-15.  Both  experimental  sequences  were  then 
repeated  at  low  (wavenumber  42)  resolution.  A  third  experiment  set  was  created  by  doubling 
the  horizontal  resolution  in  both  latitude  and  longitude,  and  repeating  the  17  forecasts  for  both 
the  ECMWF  and  NCEP-NCAR  initializations. 

The  results  are  summarized  in  Fig.  5.5,  which  is  in  the  same  format  as  Fig.  4.1a,  and 
displays  the  relative  contributions  of  various  model/data  configurations  to  the  forecast  dis¬ 
cussed  in  Chapter  4.  Each  curve  in  Fig.  4.1a  is  repeated  here,  and  three  new  curves  are  added: 
the  sensitivity  to  the  doubled  initial  uncertainty  in  waves  0-15  (open  triangles),  the  sensitivity 
to  the  doubled  initial  uncertainty  in  waves  16-42  (closed  triangles),  and  the  sensitivity  to  the 
initial  state  for  the  doubled  horizontal  resolution  (plus  signs). 

The  most  important  change  comes  from  doubling  our  prior  estimate  of  uncertainty  in 
waves  0-15.  Doubling  uncertainty  in  waves  16-42  is  also  important,  but  does  not  produce  as 
much  influence  at  longer  forecast  times  as  is  found  in  modifications  of  the  longer-wave  group. 
Doubling  the  model  resolution  also  enhances  the  error  growth  rates. 

None  of  these  experiments  produces  as  much  eight-day  (192  h)  forecast  sensitivity  to 
the  initial  state  as  the  model  selection  (xs).  As  noted  previously,  the  relative  magnitude  and 
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Fig.  5.5.  Time  evolution  of  global,  case-averaged  root-mean  square  500-mb  heights,  for  the 
1993  cases.  The  top  curve  (solid)  depicts  the  differences  of  the  analyzed  height  at  the  time  indi¬ 
cated  on  the  abscissa  minus  the  value  at  0  hours.  The  second  curve  (open  circles)  depicts  the 
differences  of  the  MRF  forecast  height  at  the  indicated  time  minus  the  analysis  at  0  hours.  The 
third  curve  (closed  circles)  depicts  the  differences  of  the  UGM  (initialized  with  NCEP)  height 
at  the  indicated  time  minus  the  analysis  at  0  hours.  The  fourth  curve  (open  squares)  depicts  the 
rms  error  of  the  height  field  forecasted  by  the  MRF.  The  fifth  curve  (closed  squares)  depicts  the 
rms  error  of  the  height  field  forecasted  by  the  UGM  (initialized  with  NCEP).  The  sixth  curve 
(xs)  depicts  the  difference  of  the  UGM  (initialized  with  NCEP)  forecast  height  minus  the  MRF 
forecast  height,  and  measures  the  sensitivity  of  the  forecasts  to  model  used.  The  seventh  curve 
(open  triangles)  depicts  the  NCEP-initialized  UGM  forecast  difference  due  to  doubled  initial 
uncertainty  in  wavenumbers  0-15.  The  eighth  curve  (plus  signs)  depicts  the  difference  of  the 
high-resolution  UGM  forecast  initialized  with  NCEP  minus  that  of  the  high-resolution  UGM 
initialized  with  ECMWF,  and  measures  the  sensitivity  of  the  high-resolution  forecasts  to  the 
initial  dataset  used.  The  ninth  curve  (closed  triangles)  depicts  the  NCEP-initialized  UGM  fore¬ 
cast  difference  due  to  doubled  initial  uncertainty  in  wavenumbers  16-42.  The  tenth  and  bottom 
curve  (diamonds)  depicts  the  difference  of  the  low-resolution  UGM  forecast  initialized  with 
NCEP  minus  that  of  the  low-resolution  UGM  initialized  with  ECMWF,  and  measures  the  sen¬ 
sitivity  of  the  low-resolution  forecasts  to  the  initial  dataset  used. 
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structure  of  the  errors  of  the  MRF  and  Utah  model  appear  to  be  similar  through  eight  days,  but 
their  forecast  differences  are  larger  than  the  sensitivity  to  presently  estimated  measures  of  the 
initial  state  uncertainty.  This  conforms  with  the  conclusion  of  Reynolds  et  al.  (1994)  for  the 
Tropics,  but  conflicts  with  their  findings  for  the  extratropics.  Sensitivity  to  initial  state  uncer¬ 
tainty  may  also  be  a  function  of  model  complexity  and  to  be  fair  in  the  relative  weighting  of 
initial  condition  differences  versus  model  differences,  the  initial  condition  sensitivity  experi¬ 
ments  should  be  done  with  a  more  complex,  higher-resolution  model. 


CHAPTER  6 


RESOLUTION  ENHANCEMENT  IN  UGM 

The  UGM  sensitivity  to  initial  state  uncertainty  was  shown  to  increase  with  resolution 
enhancement  in  the  last  chapter.  The  purposes  of  the  present  chapter  are  to  first  describe  the 
sensitivity  of  UGM  forecast  accuracy  to  both  model  and  initial  state  resolution,  and  second,  to 
repeat  at  doubled  horizontal  resolution  earlier  experiments  in  which  the  initial  state  uncertainty 
was  spectrally  binned  in  short-wave  and  long-wave  components. 

Model 

The  higher-resolution  UGM  doubles  east-west  resolution  to  include  257  grid  points. 
Fourier  series  depicting  longitudinal  variability  are  now  truncated  at  wavenumber  84  to  avoid 
aliasing  instability,  and  in  the  north-south  direction  the  number  of  model  nodal  points  is 
increased  from  82  to  164.  The  number  and  location  of  vertical  levels  is  left  unchanged  with 
respect  to  the  wavenumber  42  integrations.  Ideally,  resolution  in  all  three  dimensions  should  be 
at  least  doubled,  but  the  presently  enhanced  resolution  represents  the  largest  model  that  could 
be  executed  for  the  required  number  of  cases  within  the  timeframe  of  the  study.  Even  so,  com¬ 
puter  time  limitations  restricted  the  high  resolution  (wavenumber  84)  experiments  to  a  subset 
of  cases  reported  in  earlier  sections.  This  subset  includes  15  cases  between  26  January  and  24 


February  2004. 
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Initial  States 

A  series  of  experiments  was  performed  using  the  NCAR-NCEP  reanalyses  as  well  as 
the  GDAS  analyses  for  the  15  selected  cases.  Enstrophy  spectra  for  these  two  analyses  are  dis¬ 
played  in  Fig.  6.1  for  a  particular  day  in  the  Northern  Hemisphere  winter  2002/03.  Enstrophy 
is  the  square  of  the  vorticity  field,  and  its  spectrum  emphasizes  shorter  waves.  The  streamfunc- 
tion  is  first  computed  globally  from  the  wind  field,  and  then  projected  onto  spherical  harmon¬ 
ics.  Globally  integrated  enstrophy  may  be  expressed  as: 


(6.1) 


where  A  is  the  global  area,  T  is  the  streamfunction,  and  ^4^  is  the  amplitude  coefficient  of  the 


spherical  harmonic  component  of  degree  n  and  order  m.  The  quantity 


1/2 


is  plotted 


against  global  wavenumber  n  in  Fig.  6.1,  for  30  January  2003  at  sigma  level  0.3.  The  highest 
values  are  evident  in  the  lowest  wavenumbers,  or  in  the  longer  waves  of  both  analyses.  The 
GDAS  analysis  (circles)  clearly  retains  more  amplitude  in  global  wavenumbers  on  the  order 
of,  or  exceeding  30,  suggesting  that  it  resolves  much  more  of  the  subsynoptic  structure  than 
does  the  NCEP-NCAR  reanalysis  (solid). 


Anomaly  Correlations 

A  15 -case  subset  of  experiments  from  boreal  winter  2002/03  was  repeated  using  the 
higher-resolution  GDAS  analyses  to  initialize  the  high-resolution  (wavenumber  84)  UGM,  and 
results  are  compared  to  the  MRF  and  the  low-resolution  (wavenumber  42)  UGM  anomaly  cor¬ 
relations.  Since  this  sample  is  too  short  to  provide  an  adequate  period  for  a  priori  bias  removal. 


(1/s) 
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sqrt(Enstrophy)  03013012,  Sig  .2 


Fig.  6.1.  The  square  root  of  the  globally-averaged  enstrophy  spectrum  as  a  function  of  wave- 
number,  at  1200  UTC  on  30  Jan  2003,  for  the  NCEP-NCAR  reanalysis  (solid)  and  GDAS 

analysis  (circles)  at  sigma  level  0.2.  Units  are  s'*. 
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as  was  done  in  Chapter  3,  the  bias  is  removed  from  all  model  forecasts  using  the  average  error 
of  the  15-case  sample  (a  posteriori  bias  correction).  This  procedure  does  not  allow  a  true  mea¬ 
sure  of  forecast  accuracy,  but  may  still  be  useful  to  determine  how  much  relative  benefit  may 
be  assigned  to  improved  initial  state  analyses  compared  to  model  resolution  enhancement. 

Figure  6.2  displays  the  500-mb  height  anomaly  correlations  obtained  for  the  15-case  a 
posteriori  bias-corrected  MRF  model  forecasts  for  boreal  winter  2002/03  (squares)  and  the  a 
posteriori  bias-corrected  MRF  forecasts  obtained  for  the  17-case  sample  of  1993  MRF  predic¬ 
tions  (circles).  Anomaly  correlations  of  the  2003  MRF  model  dip  to  0.6  after  204  h,  while  they 
drop  to  this  amount  approximately  48  h  earlier  for  the  1 993  cases.  Note  that  the  curve  for  2003 
anomaly  correlations  is  different  from  that  shown  in  Fig.  3.1  due  to  the  different  case  sample 
size.  The  additional  predictive  skill  for  the  2003  MRF  model  shown  in  Fig  6.2  must  be  due  to 
either  better  initial  condition  specification  for  the  15-case  subset  relative  to  the  full  52-case 
sample,  or  the  presence  of  large-scale  atmospheric  conditions  better  suited  to  predictability 
during  the  timeframe  of  the  15-case  sample. 

As  suggested  in  Chapter  3,  the  overall  MRF  forecast  enhancement  from  1993  to  2003 
may  be  due  to  a  higher  resolution  model,  to  improved  observations,  or  to  a  possibly  more  pre¬ 
dictable  situation  in  the  2003  cases.  The  latter  explanation  may  be  unlikely,  because  the 
amount  of  skill  improvement  (48  h)  is  more  than  might  be  expected  from  the  episodic  variabil¬ 
ity  of  skill  suggested  in  Fig.  5.1  for  the  cases  studied  presently. 

The  relative  roles  of  resolution  enhancement  and  initial  state  enhancement  can  be  better 
separated  using  the  UGM.  Figure  6.3  shows  a  posteriori  bias-corrected  anomaly  correlations 
for  500-mb  geopotential  height  for  the  1993  low-resolution  UGM  runs  (initialized  with  NCEP- 
NCAR  reanalyses),  and  for  the  2003  low-  and  high-resolution  UGM  runs  (initialized  with 


Fig.  6.2.  Time  evolution  of  global,  case-averaged,  bias-corrected  anomaly  correlations  of  the 
500-mb  geopotential  height  for  the  MRF,  for  the  1993  cases  (circles)  correlated  with  NCEP- 
NCAR  reanalyses  and  2003  cases  (squares)  correlated  with  GDAS  analyses. 
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Fig.  6.3.  Time  evolution  of  global,  case-averaged,  bias-corrected  anomaly  correlations  of  the 
500-mb  geopotential  height  for  the  UGM,  for  the  1993  cases  (circles)  correlated  with  NCEP- 
NCAR  reanalyses,  and  2003  cases  at  wavenumber  42  resolution  (open  squares)  and  at  wave- 
number  84  resolution  (closed  squares),  each  correlated  with  GDAS  analyses. 
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GDAS  analyses).  The  UGM  forecasts  for  the  17-case  sample  for  1993,  executed  at  wavenum¬ 
ber  42  resolution  (circles),  produce  anomaly  correlations  exceeding  0.6  for  the  first  154  h. 
UGM  forecasts  for  the  15-case  subset  from  boreal  winter  2002/2003,  run  at  the  lower  resolu¬ 
tion  (open  squares),  produce  the  same  level  of  skill  to  approximately  170  h,  and  doubling 
model  horizontal  resolution  (closed  squares)  for  the  2003  cases  extends  the  period  of  skillful 
forecasts  by  about  another  3  h.  Again,  these  values  are  different  from  those  discussed  in  Chap¬ 
ter  3  due  to  the  difference  in  low-resolution  model  sample  size  used  for  comparison  with  the 
limited  number  of  high-resolution  runs.  Results  for  this  smaller  subset  of  cases  suggest  that 
enhancements  in  initial  conditions  from  NCEP-NCAR  reanalyses  of  1993  to  GDAS  analyses 
of  2003  play  a  substantial  role  in  predictive  skill  improvement  over  the  past  decade. 

Spectral  Binning  of  Initial  State  Uncertainty 

Initial  state  sensitivity  experiments  of  Chapter  5  are  now  repeated  using  NCEP-NCAR 
reanalyses  and  the  higher-resolution  GDAS  analyses  as  initial  conditions  for  the  wavenumber 
84  UGM.  In  view  of  the  computational  expense  of  the  wavenumber  84  model,  the  present  sec¬ 
tion  presents  experiments  in  which  the  initial  state  uncertainty  is  retained  only  in  waves  0-15, 
and  in  the  complementary  wavenumbers,  1 6-84,  rather  than  all  four  spectral  “bins”  used  for  the 
low-resolution  experiments. 

The  importance  of  higher  resolution  GDAS  initialization  to  the  wavenumber  42  UGM 
experiments  is  illustrated  in  Fig.  6.4,  which  shows  the  time  evolution  of  sigma=0.53  global 
variances  of  forecast  meridional  flow  differences  produced  by  initial  data  changes  for  the  1993 
cases  (Fig.  6.4a)  and  2002/2003  cases  (Fig.  6.4b).  For  comparison  with  2003  results.  Fig.  6.4a 
repeats  the  1993  results  displayed  in  Chapter  5,  using  initial  analysis  changes  restricted  in  the 
complementary  wave  groups  0-15  (open  circles)  and  16-42  (closed  circles),  normalized  by  the 
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Fig.  6.4.  Time  evolution  of  normalized  global  variances  at  a  =0.53,  of  forecast  meridional  flow 
differences  from  the  wavenumber  42  UGM,  produced  by  initial  data  changes  for  global  wave- 
numbers  0-15  (open  circles)  and  for  global  wavenumbers  16-42  (closed  circles).  Initial  data 
changes  in  (a)  are  from  NCEP-NCAR  to  ECMWF  reanalyses  for  the  1993  cases,  and  initial 
data  changes  in  (b)  are  from  NCEP-NCAR  to  GDAS  analyses  for  the  2003  cases. 
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effect  of  changing  the  initial  analysis  in  the  completely  resolved  spectrum.  Recall  that  the  1993 
cases  utilized  NCEP-NCAR  and  ECMWF  reanalyses.  Fig.  6.4a  demonstrates  that  the  initial 
state  meridional  flow  uncertainty  is  dominated  by  short  waves,  wavenumbers  16-42,  but  that 
the  uncertainty  of  the  longer  waves,  wavenumbers  0-15,  becomes  more  important  after  12  h, 
and  remains  more  important  for  the  14-day  (336  h)  duration  of  the  experiment.  This  is  consis¬ 
tent  with  the  findings  of  MMR 

Figure  6.4b  is  based  on  the  same  model  configuration  (UGM,  wavenumber  42),  but 
presents  results  using  the  difference  of  NCEP-NCAR  reanalyses  with  GDAS  analyses  for 
boreal  winter  2002/03  to  calculate  the  initial  state  uncertainty.  The  conclusions  are  rather  simi¬ 
lar  to  those  of  the  1993  cases  shown  in  Fig.  6.4a,  with  the  uncertainty  of  the  longer  waves 
becoming  more  important  after  12  h,  but  the  relative  importance  of  short  wave  uncertainty 
begins  to  approach  that  of  the  longer  waves  towards  the  end  of  the  14-day  forecast  period. 

Figure  6.5  shows  the  evolution  of  normalized  global  variances  at  sigma=0.53  of  fore¬ 
cast  meridional  flow  differences  produced  by  initial  data  changes  for  the  high-resolution 
research  model  for  1993  and  2003  cases.  Figure  6.5a  is  based  on  the  17-case  sample  of  1993 
experiments  that  are  identical  to  those  in  Fig.  6.4a,  but  the  results  are  based  on  the  wavenum¬ 
ber  84  UGM.  An  important  modification  of  the  results  due  to  the  increased  model  resolution  in 
Fig  6.5a  is  that  the  short-wave  uncertainty  impact  exceeds  the  long-wave  impact  for  the  first  24 
h,  compared  with  12  h  in  Fig  6.4a,  and  although  the  role  of  long-wave  uncertainty  is  most 
important  after  this  in  the  wavenumber  84  experiments,  it  is  not  as  dominant  as  in  the  wave- 
number  42  experiments.  Figure  6.5b,  which  repeats  the  2002/03  experiments  of  Fig  6.4b  but 
with  wavenumber  84  model  resolution,  suggests  that  the  influence  of  the  shorter  wave  detail 
available  in  GDAS  analyses  is  more  important  than  longer  wave  detail  for  the  first  three  days 


Fig.  6.5.  Time  evolution  of  normalized  global  variances  at  a=0.53,  of  forecast  meridional 
flow  differences  from  the  wavenumber  84  UGM,  produced  by  initial  data  changes  for  global 
wavenumbers  0-15  (open  circles)  and  for  global  wavenumbers  16-84  (closed  circles).  Initial 
data  changes  in  (a)  are  from  NCEP-NCAR  to  ECMWF  reanalyses  for  the  1993  cases,  and  ini¬ 
tial  data  changes  in  (b)  are  from  NCEP-NCAR  to  GDAS  analyses  for  the  2003  cases. 
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of  prediction,  and  the  longer  wave  uncertainty  is  more  important  after  this  time,  but  only  by  a 
relatively  small  amount. 

The  high-resolution  experiments  indicate  that  the  UGM  executed  at  wavenumber  84, 
incorporating  added  detail  available  from  the  GDAS  analyses,  displays  significantly  more 
forecast  sensitivity  to  smaller  scales  of  the  initial  state,  particularly  in  the  first  three  days  of 
prediction,  but  that  error  evolution  is  slightly  more  sensitive  to  uncertainties  in  the  larger  scales 
for  medium-  and  long-range  prediction.  These  results  modify  those  suggested  by  MMP  and  by 
experiments  described  in  Chapter  5. 

Recent  work  reported  by  Tribbia  and  Baumhefner  (2004)  also  highlights  the  need  for 
continued  research  into  the  relationship  between  uncertainty  scale  and  error  growth.  Their 
investigations  of  scale  interactions  and  short-term  predictability  with  the  NCAR  Community 
Climate  Model  Version  3  showed  that  error  growth  was  not  described  by  an  inverse  cascade, 
but  that  initial  state  uncertainty  in  the  synoptic  scales  (wavenumbers  10-20)  was  most  impor¬ 
tant  for  growth  of  errors.  Although  Tribbia  and  Baumhefner’s  experiments  used  a  highly 
sophisticated  climate  model  run  at  higher  resolutions  than  the  UGM,  qualitatively  the  results 
match  those  reported  here,  and  point  to  the  importance  of  detailed  broad  spectrum  analyses  for 


global  weather  prediction. 


CHAPTER  7 


EULER  MODEL 

Virtually  all  global  models  that  have  been  applied  for  long  periods  to  real-time  predic¬ 
tion  or  to  climate  simulations  are  based  upon  the  primitive  equations.  These  equations  filter  out 
sound  waves  and  make  the  hydrostatic  approximation.  Such  approximations  are  not  made  in 
so-called  "Euler  models"  that  are  being  developed  experimentally  in  a  number  of  limited  area 
applications  as  well  as  in  global  versions  (e.g.,  Semazzi  et  al.  1995;  Dudhia  and  Bresch  2002). 
Although  some  comparisons  have  been  made  of  results  using  the  Euler  equations  relative  to 
those  produced  by  the  primitive  equations  in  case  studies,  systematic  intercomparison  of  a 
series  of  forecasts  by  these  different  dynamical  model  cores  has  not  yet  been  publicized.  The 
present  chapter  describes  the  development  of  an  Euler  model  that  will  be  used  in  Chapter  8  to 
describe  sensitivity  to  model  dynamical  core. 

The  Euler  model  used  in  present  testing  is  relatively  new,  and  has  not  been  previously 
documented.  Equations  for  the  model  dynamics  are  presented  in  Appendix  A.  Since  the  model 
physics  is  identical  to  that  used  in  the  UGM,  and  is  documented  elsewhere  (e.g.,  Nicolini  et  al. 
1993),  its  description  is  not  repeated  here. 

Most  global  primitive  equation  models  incorporate  dynamical  cores  that  predict  vortic- 
ity,  divergence,  and  temperature,  and  diagnose  vertical  motion  and  geopotential  height  (i.e.,  the 
mass  distribution).  Vorticity  and  divergence,  rather  than  wind  components,  are  predicted 
because  the  former  variables  are  scalar  quantities  more  suited  for  expansion  in  terms  of  spheri- 
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cal  harmonics.  Such  expansions  allow  accurate  “spectral”  methods  to  be  used  in  the  approxi¬ 
mation  of  horizontal  derivatives.  Spectral  methods,  in  principle,  compute  horizontal 
derivatives  without  any  truncation  error,  aside  from  the  error  which  is  implicit  in  the  truncated 
representation  of  the  expansion.  The  UGM  applies  spectral  approximations  to  longitude  deriv¬ 
atives  and  finite  element  approximations  to  latitude  derivatives,  while  the  MRF  (GFS)  uses 
spherical  harmonic  expansion  for  all  variables  involved  in  its  dynamical  core. 

By  contrast,  most  Euler  models  use  dynamical  cores  that  predict  the  horizontal  wind 
components,  u  and  v,  as  well  as  vertical  velocity,  w,  the  temperature  or  potential  temperature, 
and  the  pressure  or  density  fields  (i.e.,  the  mass  distribution).  The  present  Euler  model  predicts 
M,  V,  w,  potential  temperature,  and  pressure,  as  outlined  in  Appendix  A.  Since  vector  fields  are 
not  well-suited  for  expansion  in  spherical  harmonics,  Euler  models  typically  use  grid  point, 
finite  difference  methods  that  retain  some  truncation  error,  even  for  resolved  structures. 

The  present  Euler  model  uses  the  “b-grid”  described  in  Appendix  A.  This  is  similar  to 
the  approach  of  Dudhia  and  Bresch  (2002),  but  differs  from  most  other  global  Euler  models, 
which  use  the  “c-grid”  (e.g.,  Semazzi  et  al.  1995).  As  pointed  out  by  Dudhia  (1993),  the  “b- 
grid”  requires  more  horizontal  averaging  to  calculate  advection  terms,  but  its  use  allows  a 
larger  time  step  than  does  the  “c-grid”  for  a  given  horizontal  grid  size. 

The  version  of  the  Euler  model  used  in  this  study  uses  a  spherical  grid  and  finite  ele¬ 
ment  approximations  for  all  flux  terms  in  the  horizontal  and  vertical  directions.  This  approach 
is,  in  theory,  more  accurate  than  most  other  second-  or  fourth-order  schemes,  and  expansion  of 
nonlinear  transports  in  terms  of  their  finite  element  projections  on  piecewise  linear  functions 
guarantees  conservation  of  certain  one-dimensional  quadratic  norms.  This  may  make  the 
method  more  stable  than  some  other  finite  difference  methods  that  do  not  impose  formal  con- 
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servation  constraints  on  second-order  moments,  and,  therefore,  do  not  contain  built-in  con¬ 
straints  to  eliminate  nonlinear  aliasing  instabilities  in  nonlinear  advection  terms. 

For  present  grid  sizes,  there  is  no  formal  reason  to  expect  substantially  different  fore¬ 
casts  in  the  Euler  model  than  in  the  UGM,  because  effects  of  vertical  acceleration  and  of  com¬ 
pressibility  are  unimportant  for  the  resolved  scales.  Some  differences  may  be  expected  in 
individual  forecasts  at  longer  lead  times  because  of  the  different  paths  along  which  chaotic 
components  of  the  flow  may  evolve,  but  there  is  little  reason  to  anticipate  systematically  dif¬ 
ferent  behaviors  between  the  UGM  and  an  Euler  model  executed  at  similar  precision.  Never¬ 
theless,  some  systematic  differences  do  emerge,  as  will  be  seen  in  subsequent  sections.  The 
following  discussion  outlines  some  inherent  model  differences  that  may  help  to  explain  differ¬ 
ing  results.  There  are  three  likely  explanations:  the  Euler  model  and  the  UGM  use  different 
vertical  resolution,  initialization,  and  turbulent  diffusion. 

The  Euler  model  requires  about  30  above-ground  vertical  levels  to  produce  stable  two- 
week  solutions,  while  the  UGM  can  be  run  with  fewer  levels,  using  only  20  above-ground  lev¬ 
els  in  present  experiments.  The  Euler  model  treats  vertical  sound  wave  propagation  implicitly, 
and  slows  down  vertical  sound  wave  propagation  in  the  process.  The  resolution  required  for 
this  process  appears  to  demand  a  more  detailed  vertical  treatment  than  needed  in  hydrostatic 
models  such  as  the  UGM,  which  filter  vertically-propagating  sound  waves  explicitly  in  the 
dynamical  formulation. 

Although  higher  vertical  resolution  places  a  greater  computational  burden  on  the  Euler 
model  than  on  the  UGM,  it  also  produces  better  resolution  of  vertical  atmospheric  variation  for 
depiction  of  vertical  wind  shear,  stratification,  and  vertical  propagation  of  resolved  phenom¬ 
ena.  This  extra  vertical  resolution  should  favor  the  Euler  model  in  forecast  accuracy,  particu- 
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larly  in  the  stratosphere,  as  outlined  below. 

It  is  also  possible  that  the  location  of  the  levels  makes  a  difference  in  forecast  outcome, 
although  this  has  not  been  demonstrated  in  the  present  study.  The  UGM  uses  a  sigma  coordi¬ 
nate  in  the  vertical,  with  levels  situated  at  sigma=0.99,  0.98,  0.96,  0.94,  0.91,  0.88,  0.84,  0.78, 
0.72,  0.66,  0.59,  0.53,  0.47,  0.41,  0.34,  0.28,  0.22,  0.16,  0.09,  and  0.03.  The  Euler  model  uses  a 
height-based  coordinate  that  carries  information  at  the  following  levels  (above  oceans):  z=0, 
100,  500,  1000,  2000,  3000,  4000,  5000,  6000,  7000,  8000,  9000,  10000,  11000,  12000, 
13000,  14000,  15000,  16000,  17000,  18000,  19000,  20000,  21000,  22000,  23000,  24000, 
26000,  28000,  and  30000  meters.  The  UGM  devotes  more  levels  (about  six)  to  the  atmospheric 
boundary  layer,  while  the  Euler  model  only  carries  three  levels  in  the  lowest  kilometer.  The 
UGM  carries  about  16  levels  in  the  troposphere,  and  only  about  four  levels  in  the  stratosphere, 
while  the  Euler  model  carries  a  similar  number  of  levels  in  the  troposphere  but  has  about  10 
more  levels  than  does  the  UGM  in  the  stratosphere. 

The  second  inherent  difference  between  the  UGM  and  Euler  model  which  may  explain 
inconsistencies  in  forecasts  from  the  two  models  is  model  initialization.  The  UGM  is  initial¬ 
ized  through  a  high-accuracy  global  interpolation  scheme  in  which  all  variables  available  from 
the  GDAS  or  reanalysis  archive  are  first  interpolated  onto  the  UGM  grid  using  spherical  har¬ 
monic  representation.  This  process  is  more  consistent  with  the  way  dynamical  variables  are 
structured  in  the  UGM,  and  the  high-order  interpolation  procedure  produces  superior  forecasts 
compared  to  simple  linear  interpolation  from  the  GDAS  or  reanalysis  grid  to  the  UGM  grid 
(preliminary  results  indicate  the  high-order  interpolation  provides  a  6  h  improvement  in  pre¬ 
dictive  skill  over  simple  linear  interpolation).  The  current  version  of  the  Euler  model  simply 
uses  low-order,  local  linear  interpolation  to  obtain  the  initial  state,  so  initialization  procedures 
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should  favor  the  UGM. 

The  final  difference  between  the  UGM  and  Euler  model  to  which  dissimilarities  in 
resulting  forecasts  may  be  attributed  is  the  treatment  of  turbulent  diffusion.  The  horizontal  dif¬ 
fusion  in  the  UGM  is  treated  analytically  in  Fourier  space  on  individual  latitude  circles.  This 
allows  implicit  time-differencing  for  diffusion,  but  makes  it  difficult  to  specify  diffusion  in  a 
way  that  depends  locally  upon  predicted  gradients.  Consequently,  the  horizontal  diffusion  is 
simply  proportional  to  the  horizontal  grid  size.  For  the  UGM  model,  the  horizontal  diffusivity 

for  vorticity,  temperature,  and  moisture  prediction  is  10^  m^  s'*,  and  diffusivity  for  divergence 
prediction  is  two  orders  of  magnitude  larger.  The  large  divergence  diffusivity  is  required 
because  of  the  time  splitting  technique  used  to  treat  gravity  waves  separately  from  advective 
and  Rossby  wave  dynamics.  The  time  splitting  provides  much  time  savings  in  the  calculation, 
and  produces  UGM  predictions  in  about  one-fifth  the  time  required  for  Euler  model  predic¬ 
tions,  but  it  typically  requires  added  diffusion  in  the  divergence  prediction. 

In  contrast,  the  Euler  model  imposes  horizontal  diffusivity  that  is  locally  proportional 
to  the  wind  field  gradients,  multiplied  by  the  square  of  the  grid  size.  The  maximum  value  is 
constrained  to  be  less  that  1.6x10^  m^  s'*  in  present  applications.  An  even  smaller  limit  is 
imposed  on  moisture  diffusion  and  thermal  diffusion.  The  results  of  the  model  are  almost  inde¬ 
pendent  of  the  value  of  diffusion  in  this  range,  signifying  that  the  Euler  solutions  are  very 
slightly  affected  by  arbitrarily  imposed  mixing  coefficients. 

The  selection  of  horizontal  diffusion  probably  favors  the  Euler  model  relative  to  the 
UGM.  The  small  values  of  diffusivity  in  the  Euler  model  allow  stronger  gradients  of  moisture 
and  stronger  moisture  convergence  in  Euler  simulations  than  in  UGM  simulations.  The  result 
is  different  precipitation  distributions,  which  are  more  realistic  in  their  areal  coverage  than  are 
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those  for  the  UGM.  Such  differences  in  precipitation  patterns  imply  differences  in  cloud  cover¬ 
age  between  the  two  models,  with  related  differences  in  the  models’  radiation  balances. 

It  is  possible  that  some  of  the  systematic  forecast  differences  that  are  found  between  the 
Euler  model  and  the  UGM  are  due  to  these  different  precipitation  and  latent  heating  rates,  par¬ 
tially  explained  by  the  different  amounts  of  diffusive  mixing.  In  view  of  the  computational 
expense  of  the  Euler  model  experiments,  it  was  not  possible  to  separate  the  effects  of  these  dif¬ 
ferent  model  characteristics  in  a  set  of  controlled  experiments.  If  the  forecast  differences  are 
indeed  produced  by  differences  in  rainfall  rates  associated  with  disparate  treatments  of  model 
diffusivities,  it  would  be  appropriate  to  characterize  them  to  be  due  to  both  dynamical  core  dif¬ 
ferences  as  well  as  physical  processes.  More  work  is  required  to  determine  the  true  source  of 
forecast  differences.  The  goal  of  the  following  results  is  merely  to  point  out  the  degree  of  fore¬ 
cast  divergence  found  in  a  systematic  one-month  intercomparison. 


CHAPTER  8 


EULER  CORRELATIONS  AND  BIAS 

Comparison  of  UGM  and  MRF  statistics  in  previous  chapters  implied  a  correlation  of 
error  patterns  and  biases  between  the  models,  despite  various  differences  in  model  configura¬ 
tion  and  resolution.  To  further  investigate  the  implication  of  a  common  limitation  to  predictive 
skill  in  forecast  models,  a  30-case  sample  of  Euler  model  forecasts  for  boreal  winter  2002/03  is 
examined  with  respect  to  results  shown  in  Chapter  3. 

Figure  8.1  shows  case-averaged,  globally-averaged,  500-mb  geopotential  height  anom¬ 
aly  correlations  for  the  various  models  without  (Fig  8.1a)  and  with  (Fig  8.1b)  an  a  posteriori 
bias  correction.  All  curves  shown  are  averages  for  the  30-cases  prior  to  26  February  2003.  Fig¬ 
ure  8.1a  indicates  that  the  Euler  model  (triangles)  has  approximately  6  h  more  predictive  skill 
than  the  UGM  (circles)  during  the  first  five  days  (120  h),  and  is  more  similar  to  the  evolution 
of  the  MRF  (closed  squares)  than  the  real  atmosphere.  Due  to  the  similarities  between  the 
UGM  and  Euler  model  (same  model  physics  and  resolution),  their  model-to-model  correlation 
is  very  high,  and  does  not  drop  below  0.6  throughout  the  14-day  (336  h)  forecast. 

Because  the  Euler  model  is  computationally  expensive,  the  completion  of  a  “training 
period”  of  cases  used  for  an  a  priori  bias  correction  was  not  possible.  Bias  corrections  are, 
therefore,  a  posteriori,  which  does  not  allow  a  true  measure  of  forecast  accuracy,  but  is  useful 
to  determine  relative  accuracies  between  the  models.  Figure  8.1b  shows  the  bias-corrected 
results  for  the  anomaly  correlations  shown  in  Fig.  8.1a.  The  Euler  model  clearly  does  not  ben- 
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Fig.  8.1.  Time  evolution  of  global,  case-averaged  anomaly  correlations  of  the  500-mb  geopo¬ 
tential  height,  for  the  2003  cases  of  the  UGM  (open  circles),  MRF  (open  squares),  and  Euler 
model  (triangles),  each  correlated  with  analyses;  and  for  the  Euler  model  correlated  with  the 
UGM  (closed  circles)  and  the  Euler  model  correlated  with  the  MRF  (closed  squares).  Correla¬ 
tions  are  (a)  not  bias  corrected  and  (b)  a  posteriori  bias  corrected. 
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efit  as  much  from  the  bias  correction  as  does  the  UGM.  Euler  model  predictive  skill  (triangles) 
now  lags  that  of  the  UGM  (circles)  by  24  h,  and  the  bias-corrected  Euler  model  is  less  able  to 
forecast  the  bias-corrected  UGM  (closed  circles)  than  in  the  uncorrected  case  (model-to-model 
anomaly  correlations  now  drop  below  0.6  at  approximately  288  h). 

This  drop  in  model-to-model  anomaly  correlations  after  bias  removal  implies  a  similar¬ 
ity  in  bias  fields  between  the  UGM  and  Euler  model.  Figure  8.2  shows  the  500-mb  meridional 
wind  bias  field  and  departures  from  zonal  mean  height  field  at  forecast  hour  336  for  the  Euler 
model,  for  the  same  30-day  period  used  for  the  UGM  and  MRF  fields  in  Chapter  3.  Once 
again,  remarkable  similarities  in  the  bias  fields  (Fig.  8.2a  and  Fig.  3.3a,b)  are  seen,  although 
there  are  some  differences,  particularly  centered  over  Japan,  the  Great  Lakes,  and  between  the 
southern  extents  of  South  America  and  Africa.  Comparison  of  the  departures  from  the  zonal 
mean  forecasted  geopotential  height  for  the  three  models  (Fig.  8.2b  and  Fig.  3.5a,b)  again 
shows  a  stronger  similarity  between  model  forecasts  than  between  any  of  the  models  and  the 
analysis  for  the  same  period  (Fig.  3.4b).  The  Euler  model,  like  the  UGM  and  MRF,  places  a 
ridge  over  the  northern  Rockies  rather  than  off  the  Pacific  coast  of  the  United  States,  and  the 
Euler  model  also  misses  the  blocking  pattern  evident  over  central  Europe  in  the  analysis. 

To  quantify  the  similarities  in  model  bias  and  investigate  the  differences  between 
Euler-to-UGM  and  Euler-to  MRF  correlations.  Fig.  8.3  displays  the  correlations  coefficients  of 
the  error  of  the  500-mb  meridional  flow  forecast  by  the  Euler  model  correlated  with  the  error 
of  the  same  variable  forecast  by  the  MRF  (Fig.  8.3a)  and  UGM  (Fig  8.3b).  Comparing  Fig  8.3a 
with  Fig.  3.6a  shows  that  the  UGM  and  Euler  model  error  fields  are  very  similar  since  correla¬ 
tions  of  both  with  the  MRF  result  in  nearly  identical  curves.  In  each  case,  correlations  of  errors 
continue  to  rise  throughout  the  14-day  (336  h)  forecasts  despite  initial  data  uncertainties  and 
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Fig.  8.2.  Euler  model  bias  and  forecast  fields  at  forecast  hour  336,  averaged  over  30  model 
runs  prior  to  26  Feb  2003.  Contours  in  (a)  depict  500-mb  meridional  wind  bias  fields  and  con¬ 
tours  in  (b)  depict  departures  from  zonal  mean  500-mb  height  field.  Contour  interval  is  2  m  s‘* 
in  (a)  and  20  m  in  (b). 


Fig.  8.3.  Time  evolution  of  case-averaged  correlation  coefficients  (not  bias  corrected)  of  the 
500-mb  meridional  wind  error  for  30  days  preceding  26  Feb  2003,  averaged  globally  (solid), 
averaged  over  the  latitude  band  90°  -20°  S  (open  circles),  and  averaged  over  the  latitude  band 
20°  -90°N  (open  squares).  Correlations  are  for  the  Euler  model  to  the  (a)  MRF  and  (b)  UGM. 
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growing  model  chaos.  As  indicated  in  Chapter  3,  this  points  to  the  growing  similarities  in 
model  bias  between  the  UGM/Euler  model  and  the  MRF. 

Figure  8.3b  reflects  correlations  of  the  Euler  model  with  the  UGM.  The  rapid  rise  in 
correlations  during  the  first  48  h  of  the  forecasts  may  be  attributed  to  uncertainties  in  the  initial 
conditions,  which  were  the  same  for  both  models.  Past  this  point,  the  correlations  tend  to  flat¬ 
ten  at  high  values,  hinting  that  the  similarity  in  bias  between  the  UGM  and  Euler  model  is  sus¬ 
tained  throughout  the  14-day  integrations. 

The  conclusions  regarding  similarity  in  model  bias  are  quantified  in  Fig.  8.4  which 
shows  the  correlation  coefficients  of  the  bias  of  the  500-mb  meridional  flow  forecast  by  the 
Euler  model  with  that  forecasted  by  the  MRF  (Fig.  8.4a)  and  the  UGM  (Fig  8.4b).  Comparison 
of  Fig  8.4a  (Euler  to  MRF  correlations)  with  Fig.  3.6b  (UGM  to  MRF  correlations)  indicates, 
again,  that  the  Euler/UGM  model  biases  continue  to  grow  in  similarity  with  the  MRF  bias 
throughout  the  forecast,  and  are  most  similar  by  day  14. 

Figure  8.4b  shows  the  Euler  model  bias  correlations  with  the  UGM  bias.  Although  the 
model  biases  are  most  similar  in  the  Southern  Hemisphere  early  in  the  forecast,  by  day  5  (120 
h)  the  model  biases  are  highly  correlated  throughout  the  globe  and  continue  to  remain  corre¬ 
lated  during  through  the  end  of  the  forecast  period. 

The  similarity  of  error  between  the  three  models  continues  to  point  to  some  fundamen¬ 
tal  limitation  of  forecast  accuracy.  Initial  state  uncertainty  is  common  to  all  models,  and  may 
explain  the  agreement  in  error  and  bias  behavior  early  in  the  integration  period,  but  it  is  still 
unclear  to  what  the  subsequently  growing  similarities  in  bias  structure  between  the  three  differ¬ 


ent  models  may  be  attributed. 
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Fig.  8.4.  Time  evolution  of  the  correlation  coefficients  of  the  500-mb  meridional  wind  bias, 
computed  over  the  30  days  preceding  26  Feb  2003,  averaged  globally  (solid),  averaged  over 
the  latitude  band  90°-20°S  (open  circles),  and  averaged  over  the  latitude  band  20°-90°N 
(open  squares).  Correlations  are  for  the  Euler  model  to  the  (a)  MRF  and  (b)  UGM 


CHAPTER  9 


EULER  ROOT-MEAN-SQUARE  AND  FORECAST  EVOLUTION 

As  indicated  in  Chapter  7,  the  forecast  evolution  of  the  UGM  and  Euler  models  is  ex¬ 
pected  to  be  similar,  and  bias  and  error  field  correlations  in  Chapter  8  showed  this  to  be  the  case. 
However,  Chapter  7  also  described  at  least  three  differences  in  Euler  and  UGM  model  formu¬ 
lation  that  may  lead  to  systematic  differences  between  the  models.  The  purpose  of  this  chapter 
is  to  further  investigate  these  differences  through  global  rms  values,  as  described  in  Chapter  4. 

Figure  9.1  displays  the  time  evolution  of  case-averaged,  global  rms  changes  of  analyzed 
(using  GDAS  analyses)  and  forecast  variables,  as  described  in  Chapter  4,  for  30  cases  from  bo¬ 
real  winter  2002/03.  Depicted  are  global  rms  changes  of  analyzed  height  (Fig.  9.1a)  and  wind 
(Fig.  9.1b)  (solid)  as  well  as  the  global  rms  magnitude  of  the  forecast  change  for  these  variables 
for  the  MRF  (open  circles),  UGM  (open  squares),  and  Euler  model  (open  triangles).  Addition¬ 
ally,  Fig.  9.1  shows  the  forecast  error  for  the  height  and  wind  for  the  MRF  (closed  circles),  UGM 
(closed  squares)  and  Euler  model  (closed  triangles),  as  well  as  two  measures  of  sensitivity  of 
height  and  wind  forecasts  to  model  selection.  Model  sensitivity  is  measured  for  a  change  from 
the  Euler  model  to  the  UGM  (diamonds)  and  from  the  Euler  model  to  the  MRF  (xs). 

As  described  in  Chapter  4,  the  top  (solid)  curves  of  Fig.  9.1a  and  Fig.  9.1b  depict  an  ap¬ 
proximate  measure  of  the  observed  change  in  the  signal  associated  with  each  field  (i.e.,  the  error 
of  a  persistence  forecast).  Below  these  curves  are  the  rms  evolution  of  each  field  from  initial  to 
forecast  time  (model  signals)  for  each  of  the  three  models.  For  both  the  height  and  the  wind 
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Fig.  9.1.  Time  evolution  of  global,  case-averaged  rms  500-mb  heights  and  winds  for  the  boreal 
winter  2002/03  30-case  sample.  The  top  curve  (solid)  depicts  the  differences  of  the  analyzed 
variable  at  the  time  indicated  along  the  abscissa  minus  the  value  at  0  h.  The  second  curve  (open 
circles)  depicts  the  differences  of  the  MRF  forecast  variable  at  the  time  indicated  minus  the 
analysis  at  0  h.  The  third  curve  (open  triangles)  depicts  the  differences  of  the  Euler  model  fore¬ 
cast  variable  at  the  time  indicated  minus  the  analysis  at  0  h.  The  fourth  curve  (open  squares) 
depicts  the  differences  of  the  UGM  forecast  variable  at  the  time  indicated  minus  the  analysis  at 
0  h.  The  fifth  curve  (closed  triangles)  depicts  the  rms  error  of  the  variable  forecast  by  the  Euler 
model.  The  sixth  curve  (closed  squares)  depicts  the  rms  error  of  the  variable  forecast  by  the 
UGM.  The  seventh  curve  (xs)  depicts  the  difference  of  the  Euler  forecast  variable  minus  the 
MRF  forecast  variable,  and  measures  the  sensitivity  of  the  forecasts  to  the  model  used.  The 
eighth  curve  (closed  circles)  depicts  the  rms  error  of  the  variable  forecast  by  the  MRF.  The 
ninth  curve  (diamonds)  depicts  the  difference  of  the  Euler  forecast  variable  minus  the  UGM 
forecast  variable,  and  also  measures  the  sensitivity  of  the  forecasts  to  the  model  used.  Curves 
in  (a)  are  for  500-mb  heights,  and  curves  in  (b)  are  for  500-mb  meridional  wind. 


meters  per  second 


85 


field,  the  MRF  (open  circles)  maintains  slightly  less  variability  than  the  atmosphere,  except  for 
the  height  field  at  long  lead  times  (after  288  h).  The  Euler  model  (open  triangles)  and  UGM 
(open  squares)  both  retain  even  less  variability  than  does  the  MRP,  although  the  Euler  model 
signal  is  slightly  stronger  than  that  of  the  UGM  for  both  variables  at  all  forecast  times.  This  is 
consistent  with  the  weaker  diffusion  used  in  the  Euler  model  compared  to  the  UGM,  described 
in  Chapter  7  as  one  of  the  inherent  differences  between  the  two  models. 

Below  the  model  signal  curves  are  the  rms  errors  for  the  height  and  wind  for  the  MRF 
(closed  squares),  UGM  (closed  circles)  and  Euler  model  (closed  triangles).  As  seen  in  Chapter 
4,  these  curves  asymptote  toward  the  respective  model  forecasted  signals  after  day  eight  for  the 
winds,  with  slower  relative  error  growth  for  the  height  field. 

The  final  two  curves  of  Fig.  9.1  measure  the  sensitivity  of  the  height  and  wind  forecasts 
to  the  model  used.  For  the  30-case  sample,  this  measure  of  sensitivity  for  a  change  from  the  Eu¬ 
ler  model  to  the  UGM  (diamonds)  is  less  than  the  forecast  error  of  the  three  models.  This  sup¬ 
ports  the  idea  that  the  Euler  model  and  UGM,  started  from  the  same  initial  conditions,  resemble 
each  other  more  than  either  resembles  reality.  The  sensitivity  for  a  change  from  the  Euler  model 
to  the  MRF  (xs)  surpasses  (for  the  height)  or  approaches  (for  the  wind)  the  forecast  error  curves 
for  the  three  models.  Therefore,  the  Euler  model’s  ability  to  anticipate  the  evolution  of  a  more 
sophisticated,  higher-resolution  model  is  approximately  as  good  as  the  ability  of  each  of  the 
models  to  predict  the  atmosphere.  The  fact  that  the  sensitivity  for  a  switch  from  the  Euler  model 
to  the  UGM  is  nearly  two-thirds  the  value  of  the  sensitivity  for  a  switch  from  the  Euler  model 
to  the  MRF  is  evidence  that,  although  the  Euler  model  and  UGM  have  similar  bias  and  error 
fields,  there  are  systematic  differences  between  the  two  models  that  results  in  disparate  forecast 
evolution.  Determining  the  relative  contributions  of  the  three  inherent  differences  in  Euler  mod¬ 
el  and  UGM  formulation  detailed  in  Chapter  7  is  beyond  the  scope  of  the  current  study. 


CHAPTER  10 


SUMMARY  AND  CONCLUSIONS 

This  study  has  investigated  predictability  using  four  different  global  forecast  models, 
and  three  different  estimates  of  the  initial  state  selected  from  reanalyses  for  1 993  and  opera¬ 
tional  data  for  boreal  winter  2002/03.  Our  approach  has  been  to  compare  the  different  model 
forecasts  and  errors,  and  the  sensitivity  of  forecasts  to  the  model,  to  model  resolution,  and  to 
the  initial  state.  One  of  the  models  (the  experimental  UGM)  was  applied  with  all  three  analyses 
(ECMWF  and  NCEP-NCAR  reanalyses  for  1993/2003  and  GDAS  analyses  for  2003).  The 
Euler  model  was  initialized  only  with  the  GDAS  analyses  for  2003,  and  the  other  models  are 
different  versions  of  formerly  and  currently  operational  MRF  models  initialized  with  NCEP 
analyses. 

The  comparison  of  results  of  experimental  models  with  operational  models  allows 
more  predictability  inferences  than  would  the  use  of  a  single  operational  model.  The  strategy 
has  been  to  maintain  substantially  lower  resolution  and  a  simplified  physical  parameterization 
in  the  UGM  and  Euler  model  than  is  available  in  the  MRF.  This  allows  for  the  assessment  of 
the  roles  of  resolution  and  of  model  sophistication  relative  to  the  improvements  of  the  initial 
state  specification.  The  relative  contributions  of  these  separate  improvements  cannot  be  deter¬ 
mined  from  a  single  forecast  system  in  which  model  initialization  as  well  as  sophistication  are 
changing  simultaneously.  The  first  guess  is  such  an  important  component  of  the  data  assimila¬ 
tion  cycle  that  model  advances  will  lead  to  initial  state  advances.  Consequently,  the  value  of 
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individual  contributions  is  not  entirely  clear  even  in  the  present  approach,  but  the  strategy  does 
help  sharpen  the  questions  and  diagnostic  methods. 

Earlier  low-resolution  experiments  by  MMP  were  expanded  from  5  to  14  days.  Addi¬ 
tionally,  some  experiments  were  performed  with  doubled  resolution,  and  modified  estimates  of 
initial  state  uncertainty.  The  relative  magnitude  of  the  forecast  uncertainty  associated  with 
model  changes  was  compared  to  the  role  of  the  initial  state  uncertainty  estimates.  The  main 
results  are  summarized  below. 

The  results  of  the  1993  reanalysis  cases  suggest  that  model  bias  extraction  produces 
more  benefit  to  14  days  than  do  other  tested  modifications  of  the  model,  including  changes  to 
the  model  heating.  Bias-corrected  anomaly  correlations  produced  by  the  UGM  lag  MRP  bias- 
corrected  anomaly  correlations  by  about  12  h  to  day  6  (144  h),  and  become  very  similar  to  the 
MRP  bias-corrected  anomaly  correlations  by  day  7  (168  h).  Bias  corrected  anomaly  correla¬ 
tions  fall  below  0.6  after  6.5  days  for  the  UGM  initialized  with  the  NCEP-NCAR  reanalyses 
and  fall  below  0.6  after  7  days  for  the  UGM  initialized  with  ECMWF  reanalyses.  These  results 
suggest  that  the  bias-corrected  forecasts  of  a  relatively  simple  research  model  such  as  the 
UGM  have  almost  as  much  value  in  predicting  large-scale  features  as  those  obtained  from  a 
more  highly  developed  model  such  as  the  MRP  that  was  used  operationally  until  approxi¬ 
mately  1995. 

The  results  for  the  2003  real-time  sample  suggest  that  the  operational  MRP  has  gained 
approximately  36-48  h  of  forecast  value  (depending  on  sample  size)  since  the  earlier  experi¬ 
ments,  and  exhibits  bias-corrected  anomaly  correlations  of  approximately  0.6  past  8  days  for 
the  presently  studied  samples.  Meanwhile  the  UGM  executed  with  the  more  recent  data  has 
increased  the  period  of  useful  predictability  by  only  approximately  6-16  h  (depending  on  sam- 
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pie  size)  relative  to  the  1993  data,  and  gains  only  3  h  predictive  skill  when  executed  at  doubled 
horizontal  resolution.  It  is  possible  that  the  extra  30-h  gain  made  by  the  MRF  is  due  to  the 
model  advances  at  NCEP  rather  than  improvements  of  the  initial  state.  Although  this  possibil¬ 
ity  cannot  be  discounted,  the  error  patterns  of  the  low-resolution  UGM,  low-resolution  Euler 
model,  and  high-resolution  MRP  show  strong  positive  correlation  through  the  prediction,  sug¬ 
gesting  that  much  of  the  remaining  error  may  be  explained  by  some  common  deficiency  of  all 
three  models. 

The  MRP  contains  certain  inherent  advantages  in  the  balance  of  initial  fields,  since  the 
same  model  is  incorporated  to  provide  the  initial  guess  fields  for  the  analyses.  The  UGM  inter¬ 
polates  the  reanalyses  values  to  UGM  grid  points  using  spherical  harmonic  representation,  and 
the  Euler  model  uses  local  linear  interpolation.  Since  the  short-term  evolution  includes  a 
spinup  time  for  both  models  to  adjust  to  small  local  imbalances,  initialization  procedures  will 
favor  the  more  sophisticated  MRP  model.  In  view  of  these  simplifications  for  the  UGM  and 
Euler  models  it  is  somewhat  surprising  that  the  length  of  skillful  forecasts  does  not  differ  by 
more  than  roughly  12  h  for  the  1993  sample.  This  result  is  consistent  with  the  possibility  that 
those  predictions  were  not  strongly  limited  by  differences  in  model  sophistication,  but  were 
more  limited  by  uncertainty  in  the  initial  state. 

To  better  understand  the  sensitivity  to  model  heating,  all  UGM  experiments  presently 
shown  were  repeated  by  substituting  climatological  values  obtained  from  reanalysis  archives 
for  model-generated  values  of  radiative  heating  of  the  atmosphere.  All  statistics  matching 
those  displayed  in  this  paper,  as  well  as  individual  forecasts,  were  remarkably  similar  to  the 
UGM  forecasts  in  which  internally  consistent  physical  parameterizations  were  computed 
rather  than  prespecified,  and  are  not  shown.  Similar  experiments  were  made  with  changes  to 
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the  UGM  values  of  diffusion  for  divergence  and  moisture,  will  little  resulting  differences  in 
large-scale  flow  predictions. 

The  spectrally  binned  experiments  performed  by  MMP,  and  extended  here  to  14  days, 
were  also  performed  with  the  UGM  at  doubled  horizontal  resolution  using  both  1993  and  2003 
data.  These  experiments  resulted  in  modifications  of  previous  hypotheses  regarding  the  impor¬ 
tance  to  forecast  evolution  of  the  uncertainty  of  long  waves  relative  to  the  uncertainty  of  short 
waves  in  the  initial  state.  Doubled  model  resolution  incorporated  more  detail  from  the  GDAS 
analyses  at  shorter  wavelengths,  allowing  significantly  more  forecast  sensitivity  to  the  smaller 
scales  of  the  initial  state,  particularly  in  the  first  few  days  of  prediction.  However,  uncertainty 
in  the  larger  scales  of  the  initial  state  still  plays  an  important  role  in  medium-  and  long-range 
forecasts.  Such  results  point  to  the  importance  of  detailed  analyses  at  all  scales  for  improved 
global  weather  prediction. 

Our  experiments  indirectly  suggest  that  a  major  benefit  enjoyed  by  the  current  version 
of  the  MRF  is  the  extra  resolution  of  that  model  (T254  and  T170)  relative  to  the  other  models 
tested  in  the  present  experiments  (highest  resolution  of  T84).  As  previously  discussed,  higher 
resolution  also  produces  more  sensitivity  to  the  initial  state,  possibly  allowing  for  a  more  com¬ 
plete  accounting  for  the  many  initial  state  improvements  that  have  been  incorporated  into  the 
analyses  since  1993. 

The  importance  of  model  physics  to  predictive  skill  was  not  explicitly  tested  in  the 
present  study.  The  UGM  and  Euler  models  have  very  simple  physical  parameterizations  com¬ 
pared  to  those  used  in  both  versions  of  the  MRF  utilized  here.  Still,  despite  the  fact  that  Euler 
model’s  large-scale  predictive  skill  was  comparable  to  (not  bias  corrected)  or  less  than  (bias 
corrected)  that  of  the  UGM,  the  Euler  model  appears  to  evolve  more  realistic  smaller-scale 
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physical  processes,  such  as  precipitation.  Figure  10.1  is  a  comparison  of  the  average  of  14  con¬ 
secutive  days  from  boreal  winter  2003  for  the  accumulated  precipitation  at  day  11  for  the 
UGM  (Fig.  10.1a),  Euler  model  (Fig.  10.1b),  and  the  MRF  (Fig.  10.1c).  UGM  precipitation  is 
unrealistically  tied  to  the  orography  (i.e.,  the  Andes)  as  a  result  of  frictional  effects  and  the  use 
of  sigma  coordinates,  and  although  the  UGM  rains  more  than  the  Euler  model,  the  precipita¬ 
tion  from  the  UGM  tends  to  be  concentrated  in  “bulls-eyes”  throughout  the  globe.  The  areal 
coverage  of  the  Euler  model  precipitation  forecast  is  more  realistic  and  more  closely 
approaches  that  of  the  MRF.  Still,  neither  the  UGM  nor  the  Euler  model  is  close  to  producing 
the  type  of  precipitation  features  evident  in  the  more  sophisticated  MRF  model,  and  none  of 
the  models  accurately  reflects  observed  precipitation  values  for  the  period  of  interest,  shown 
for  South  America  in  Fig.  10.2. 

Although  the  models  do  exhibit  forecast  differences  and  have  varying  complexities, 
resolutions  and  dynamical  core,  all  three  produced  striking  similarities  in  500-mb  meridional 
wind  biases  by  day  14.  These  similarities  hint  at  some  common  deficiency  in  these  models 
which,  if  present  in  all  numerical  models,  would  have  important  implications  for  ensemble 
forecasting.  If  all  models  tended  to  diverge  from  reality  in  similar  ways,  ensemble  methods 
may  not  add  value  to  the  long-range  forecast  process. 

All  speculations  could  be  checked  by  repeating  the  experiments  displayed  here  using  a 
much  higher  resolution,  more  sophisticated  model  for  experimental  forecasts  and  systemati¬ 
cally  modifying  its  input,  resolution,  and  model  parameterizations  to  quantify  the  role  of  model 
resolution  and  sophistication  on  predictive  skill.  Comparison  of  the  forecast,  bias  and  error 
fields  from  a  fourth,  independent  model  with  those  from  the  UGM,  Euler  and  MRF  models 
would  also  allow  further  assessment  of  the  role  of  improvements  in  initial  state  specification. 
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Fig.  10.1.  Total  accumulated  precipitation  at  forecast  day  11,  averaged  over  14  consecutive 
days  from  28  Jan  2003  through  10  Feb  2003.  Totals  are  shown  for  the  (a)  UGM,  (b)  Euler 
model,  and  (c)  MRP  models.  Units  are  in  centimeters. 


APPENDIX  A 


EULER  MODEL  DESCRIPTION 

The  governing  equations  for  zonal,  meridional,  and  vertical  flows  (u,  v,  and  w)  and  for 
pressure,  potential  temperature,  and  specific  humidity  (p,  0 ,  and  q)  are 
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Here  X, ,  (p ,  and  r  are  longitude,  latitude,  and  the  earth’s  radius,  V  is  the  velocity  vector,/is  the 
Coriolis  parameter,  Q  is  the  angular  speed  of  the  earth’s  rotation,  p  is  the  density  of  air,  s  is 
the  material  surface  top  of  the  model,  s  is  the  initial  height  of  the  material  surface,  z*  is  the 
terrain-following  vertical  coordinate,  Zq  is  the  ground  elevation,  g  is  gravity,  Q  is  the  diabatic 

heating,  E  is  evaporation,  P  is  precipitation,  and  (where  X  is  u,  v,  w,  0 ,  or  q)  is  turbulent 

mixing.  The  terrain  slope  angle  is  neglected  for  radiation  calculations  in  present  applications. 
The  parameters  a  and  6  are  artificial  parameters  used  to  slow  down  horizontal  and  vertical 
sound  wave  propagation,  and  are  set  equal  to  1.0.  The  ratio  of  specific  heats,  y,  is  1.4. 
Turbulent  mixing  is  defined  as 


(a.7) 


where  the  horizontal  exchange  coefficient,  is  equal  to  0.36(r  x  A(p)^//(V  •  KxV)  ,  and 

is  bounded  by  (2i'A(p)^/10^  «  2.4  x  10**  m^  s‘^  The  vertical  exchange  coefficient,  K,  is  given 
by  a  level  2.5  prognostic  turbulence  treatment,  as  in  Nicolini  et  al.  (1993). 

The  diabatic  heating,  Q,  includes  radiative  and  latent  heating,  as  described  by  Nicolini 
et  al.  (1993),  and  the  thermodynamic  equation  includes  surface  heat  balance  and  subsurface 
temperature  prediction  on  five  levels  below  the  soil,  as  described  by  Paegle  and  McLawhom 
(1983). 


Advection  is  computed  as 
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The  vertical  coordinate,  z* ,  follows  the  topography  and  is  defined  by 


7*  z= 


(s-Zq) 


(A.9) 


as  in  Mahrer  and  Pielke  (1977)  (see  Figure  A.l).  If  s  equals  5=altitude  of  the  top  of  the  model 
atmosphere  (constant),  then 

w*  =  —^w-^^^iV-VZf^)  .  (A.IO) 

S-Zq  S-Zq 

Additionally,  the  pressure  perturbation,  p' ,  is  introduced  in  order  to  minimize  truncation  error 
in  the  vicinity  of  steep  terrain  in  the  calculation  of  pressure  gradient  forces.  The  total  pressure 
is  p  =  p^  +  p' ,  and  the  reference  state  pressure,  depends  only  on  z. 

The  dependent  variables  are  vertically  staggered,  as  in  Figure  A.2.  The  equations  gov¬ 
erning  p'  and  w  are  coupled  into  a  single  equation  for  the  forecasted  w  at  time  step  n+1.  The 

time  and  vertically  discretized  equation  for  p'l*^  is 


+  1  *n-\ 

Pk  '^Pk  ■>  ^ 


Pii*r(i±D  »+ 1  +  (1 
2  2  2 


Wk_i  + 


+ 


Wi 


r 


-Wk 


1 


rcos(p 


(dul  dulfdZ(p.fz* 


+ 


^dv"k  ^  dvl  rdz 


+  ( -J-J\  r  ( Q 

V5(p  /  \S  —  Z(J J  -  Z(/ 


dz*^d(py\s-z 


-W; 

2  * 


(A.11) 


where  ()*=[(  )(t+(  )*+i]/2  and  P=1  is  used  to  produce  a  fully  implicit  backward 
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Fig.  A.l.  Diagram  of  Euler  model  vertical  coordinate. 
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Fig.  A.l.  Vertical  staggering  of  dependent  variables. 
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scheme.  A  similar  equation  is  used  to  obtain  and  these  equations  are  solved  for  ' 


and  p'1^ }  which  are  inserted  in  the  vertical  momentum  equation: 
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where  p'  =  {p  -  p^RT)/{RT;) ,  T'  =  +  0.61^ r,  and  a =1.  Basic  state  ther- 

Po 

modynamic  quantities  (7^,  ,  and  p^)  are  computed  from  global  averages  of  initial  values  on 

constant  height  surfaces,  and  are  interpolated  to  terrain  following,  z* ,  coordinates.  The  final 
equation  for  w”  ^  '  is  of  the  form: 


.  «+i,r>  ''■'■•ly-'  o'’ 

^k^k-l'^^k^k  '^^k^k+\  ~  ^k 


(A.13) 


where 


(A.14) 


Bk  =  B\  + 


■{Kk-i+Kk)  ^  (Kk  +  Kk^Q- 


(A.15) 


(Ar^  + A!’^+ i)/[(z^+ j -Zyj.)(z^+ j  —  j)] 


(A.16) 


n-\ 

W.  ^  n 


sa 


"1  ““  (3 ,  1  i«“ 

)j- 


M 

ypk 


\-^ipk  ^Pk+l) 


2  2 
Ko-/pk)  (1  +  p) 

S-Zq  2(Z^+j-Z;j,) 


9^  k 


(A.I7) 


101 


n-\  ^  «  fV-[F//(5 -%)]!" 


{s-Zq) 


n-\  n-\ 


(A.18) 


A\  =  XQ\,{-XQ2,-XQ^,)+XQA,{XQ2,+XQ3,)  (A.19) 
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and  the  following  definitions  apply: 
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Note  that  XQljt,  XQ2jt,  XQ3^,  XQ4^,  and  XQ6/.  vary  horizontally  on  terrain-following  sur¬ 
faces,  as  well  as  vertically. 

The  solution  for  *  is  obtained  by  Gauss  elimination,  subject  to  boundary  condi¬ 
tions  ^=0  at  the  top  and  bottom  of  the  model.  This  filters  horizontally-propagating,  external 
gravity  waves.  The  solution  for  w  is  only  slightly  affected  by  horizontal  and  vertical  diffusion. 
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and  the  diffusion  terms  (involving  K  and  are  generally  ignored  in  present  applications. 
After  is  obtained,  these  values  are  used  in  the  equation  for  pressure  to  obtain  p'Jt  ^  •  The 
parameters  a  and  6  were  introduced  in  early  testing  of  the  model.  If  a  <  1 ,  then  vertical 
sound  wave  propagation  is  slowed  by  a  factor  Ja .  If  5  <  1 ,  then  three-dimensional  sound 


wave  propagation  is  slowed  by  the  factor  Vb  •  These  tactics  are  not  employed  in  any  of  our 
integrations. 

The  equations  for  horizontal  momentum  and  for  specific  humidity  use  centered  time 
differencing  to  approximate  d()ldt 


dt 


2At 
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and  all  terms  in  these  equations,  with  the  exception  of  horizontal  and  vertical  diffusion  are 
specified  at  time  step  n.  Horizontal  diffusion  is  specified  at  time  step  n-1  and  vertical  diffusion 
at  time  step  n+ 1. 

The  vertical  structure  of  solutions  for  u,  v,  and  q  are  expanded  in  finite  element  repre¬ 
sentation,  and  the  resulting  implicit  equations  for  u,  v,  and  q  are  of  the  form 

k+ \  ~  ^k  (A.28) 

where  ryj.=M^,  v^,  or  and  solutions  for  are  obtained  subject  to  lower  and  upper  boundary 
conditions  on  F,  and  F^^^  •  ^or  wind  components,  Fj  =0,  while  for  specific  humidity  the  lower 
boundary  condition  specifies  either  qj  at  a  water  surface  or  (dq/dz) ;  at  a  land  surface.  In  the  lat¬ 
ter  case,  the  vertical  gradient  of  q  depends  upon  surface  evaporation.  The  top  boundary  condi¬ 
tion  maintains  F^^  equal  to  the  initial  value. 
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The  equations  for  atmospheric  potential  temperature  (0)  and  for  soil  temperature  (T^) 

also  use  centered  time  differencing  and  all  terms  of  these  equations,  with  the  exception  of  hor¬ 
izontal  and  vertical  diffusion,  are  specified  at  time  step  n.  Horizontal  diffusion  is  specified  at 
time  step  n-1  and  vertical  diffusion  at  time  step  n+ 1. 

Figure  A. 3  shows  the  vertical  distribution  of  temperature  variables  and  governing 
equations  in  the  atmosphere,  interface  and  subsurface.  Here  p ,  c,  and  K  represent  density,  spe¬ 
cific  heat,  and  vertical  diffusivity  for  the  atmosphere  (subscript  yf)  or  the  soil  (subscript  s).  is 
soil  temperature  and  and  LHF  are  infrared,  solar  and  latent  heat  fluxes  at  the  surface. 

The  vertical  structure  of  0  in  the  atmosphere  and  of  7^  in  the  soil  are  expanded  in  finite 
elements,  leading  to  tri-diagonal  algebra.  The  interface  condition  uses 

Q2-Q1  _  8B 

z\-Zq  dz 

and  internally  computed  surface  infrared  flux,  F^,and  surface  solar  radiative  flux,  .  LHF  is 
currently  prespecified.  Here  K,  K^,  and  are  conductivities  of  air  (K  and  Kj)  and  soil  or 
water  (K^),  and  may  differ  from  the  diffusivity  used  in  the  other  equations.  The  surface  type  is 
specified  by  ,  c^,  and  K^,  and  varies  with  soil  type  or  water.  Roughness  length,  Zq,  is  prespec¬ 
ified,  as  is  the  subsurface  temperature  at 

The  vertical  structure  solution  for  0”^^  and  consists  of  a  single  tri-diagonal 

matrix  inversion,  extending  from  ^'kz- 

The  horizontal  variation  is  carried  on  a  “b-grid,”  with  horizontal  flow  variables 
depicted  on  a  staggered  grid  with  respect  to  other  variables  (see  Figure  A.4).  The  model  uses  a 
spherical  coordinate,  and  the  poles  carry  w,p\  0 ,  and  q. 


xfr. 


^s,  0  ~  ^s,  -1 
Zq  - 


3 

dz 


(A.29) 
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Fig.  A.3.  Vertical  distribution  of  temperature  variables  and  governing  equations  in  the  atmo¬ 
sphere,  interface  and  subsurface. 
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All  advection  terms  are  rewritten 

V-VT  =Vh-VhT  +  w^  =  Vfj- ( FF)  -  r(V^  •  Vh)  +  (A.30) 


where  F=(m,v)  and 


V//() 


1  a()>.,ia(H 

rcos(p9A.  rd(p 


(A.31) 


Calculation  of  fluxes  on  the  staggered  “b-grid”  is  done  in  the  usual  manner  by  averaging  wind 
components  and  the  fluxed  quantity  to  common  intermediate  points. 

The  model  retains  options  to  use  finite  difference  or  finite  element  (compact  implicit) 
approximations  for  flux  divergence.  The  finite  element  approach  appears  to  be  more  accurate 
and  more  stable,  and  is  used  presently. 

Convergence  of  meridians  produces  severe  time  step  restrictions.  A  fourier  filter  is 
applied  in  longitude  to  alleviate  this  problem.  The  filter  removes  all  longitudinal  waves  with 
wavelength  less  than  the  shortest  wavelength  locally  resolved  in  latitude. 

The  model  also  has  options  for  variable  latitude  grid  spacing  and  for  a  mathematical 
coordinate  rotation  that  allow  a  region  of  telescoping,  refined  resolution  within  a  selected  sub- 
domain,  similar  to  the  method  described  for  the  UGM  (primitive  equation  model)  in  Paegle 
(1989.) 

The  model  is  similar  to  that  described  by  Dudhia  and  Bresch  (2002),  but  uses  a  single, 
unified  coordinate  (rather  than  the  patched  map  projection  used  by  Dudhia  and  Bresch)  and 
allows  variable  resolution.  The  current  model  does  not  incorporate  time  splitting  for  acoustic 
modes  (as  used  by  Dudhia  and  Bresch).  It,  therefore,  requires  a  small  time  step  and  is  relatively 


computationally  slow. 


APPENDIX  B 


ASYMPTOTIC  LIMITS  FOR  CORRELATION  COEFFICIENT 

The  correlation  coefficients  of  the  error  of  the  500-mb  meridional  flow  forecast  by  the 
UGM  correlated  with  the  error  of  the  same  variable  forecast  by  the  MRF  is  calculated  by 

corrcoefF=  - s - ;  1/2  (B.l) 

where  Vu,  and  are  the  UGM-predicted,  MRF-predicted  and  analyzed  500-mb  meridi¬ 
onal  winds,  respectively,  and  (  )  represents  the  global  average.  This  measure  of  error 
correlation  asymptotes  to  a  nonzero  value  at  long  times  when  predictability  is  lost  due  to  initial 
state  uncertainty  and  to  the  chaotic  nature  of  the  atmosphere.  The  limiting  value  appears  to  be 
approximately  0.5  in  long  simulations,  and  a  similar  value  can  also  be  expected  from  the  fol¬ 
lowing  analysis. 

The  correlation  can  be  expanded  as 

C=  "-J  -1 - -  -T  ^  _  1/2  ■  (B.2) 

After  a  sufficiently  long  time,  random  errors  of  the  initial  state  would  amplify  differently  in  the 
two  models,  which  might  then  diverge  sufficiently  from  reality  and  from  each  other  so  that  the 
average  of  many  forecasts  would  imply 


(B.3) 
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In  this  case,  C  approaches  0.5  in  the  long  time  limit  average  of  many  forecasts.  Since  the  mod¬ 
els  are  not  nearly  perfect,  and  the  UGM  systematically  underestimates  the  wave  activity,  the 
precise  value  to  which  C  should  asymptote  in  the  long  time  limit  of  vanishing  predictability  is 
unclear,  but  empirical  evidence  points  to  a  value  close  to  0.5. 
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